PROGRAM iFIT C Iris Breda (email: Iris.Breda@astro.up.pt) C Affiliation: C i) Instituto de Astrofísica e Ciências do Espaço - Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal C (www.iastro.pt), and C ii) Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal C C ========================================================== C Purpose: This program fits the Sersic Law to an observed C Surface Brightness Profile using the C fitting concept described in Breda et al., 2019 (A&A, in press) C C Input parameters: C 1.filename - ex: NGC0810_SBP_g.dat C x(i) --- array with R values (arcsec) C y(i) --- array with surface brightness (SB) values (mag/arcsec2) C sig_y(i) --- array with uncertainties in SB (optional) C 2.ynPSF: (y/n) depending on whether the user wishes to correct for PSF convolution effects C . typePSF (if ynPSF = y): type of the PSF - (M) for Moffat, C (G) for Gaussian or (U) for user input PSF; C . FWHM (if typePSF = M/G): desired FWHM for the PSF; C . PSFfilename (if typePSF = U): name of the file where is the C user input PSF C x(i) --- array with R values (arcsec) C y(i) --- array with SB values (mag/arcsec2) C 3.distance: distance to the galaxy in Mpc *(or 0 instead, not C computing distance-dependent quantities) C C To Run: iFIT_FV filename ynPSF typePSF FWHM/PSFfilename distance C C ex: iFIT_FV NGC0810_SBP_g.dat y U NGC0810_PSF_1.dat 102 C ex: iFIT_FV SBP.dat n 10 C ex: iFIT_FV SBP.dat y M 1.5 10 C ex: iFIT_FV SBP.dat y U PSF.dat 0 C C Output files: C 1.filename.Results.dat - ex: NGC0810_SBP_g.dat.Results.dat C eta --- Sersic index C Reff --- effective radius (arcsec) C Reff* --- effective radius (kpc) C SBReff --- surface brightness at Reff (mag/arcsec2) C AppM --- apparent magnitude (mag) C AbsM* --- absolute magnitude (mag) C SB(R=0) --- central surface brightness (mag/arcsec2) C 2.filename.SProf.dat - ex: NGC0810_SBP_g.dat.SProf.dat C x(j) --- array with modeled R values (arcsec) C y(j) --- array with modeled SB values (mag/arcsec2) C if convolved C xC(j) --- array with modeled conv. R values (arcsec) C yC(j) --- array with modeled conv. SB values (mag/arcsec2) C 3.filename.ps - ex: NGC0810_SBP_g.dat.ps C - visual output of best fit to the data C C To compile: C gfortran -o iFIT_FV iFIT_FV.f -L/usr/local/pgplot -L/usr/X11/lib -lpgplot -lX11 C =========================================================== INTEGER m,c,i,j,g,f,k,nmax,lim(3),etamax,jstep,nlin0 INTEGER n,npoints,ier,ncol INTEGER opt_flag,INFO,LWA,IWA(2) INTEGER scycle,etalim(2),run PARAMETER (nmax=25000) CHARACTER filename*40,filenamePSF*40,yPSF*1,MGU*1,fw*10 CHARACTER string*100 REAL xstep(2),minim(250),FWHM,distance,xmax,ymin,ymax REAL eta,etamin,etaT(250),calC,chi_flag,eeta DOUBLE PRECISION xi(nmax),yi(nmax),xf(nmax),yf(nmax),ysym(nmax) DOUBLE PRECISION xmod(nmax),ymod(nmax),PSF(nmax),conv(nmax) DOUBLE PRECISION oy(nmax),xsym(nmax),ysymc(nmax) DOUBLE PRECISION Rx(10),Ix(10),lum,mag,rmin,rmax,magAbs DOUBLE PRECISION bn,val(40),slo,eslo DOUBLE PRECISION oR40,oI40,oR60,oI60,oR70,oI70,oR80,oI80 DOUBLE PRECISION R50,oR50,oI50,oR90,oI90,Mof(nmax) DOUBLE PRECISION b,eb,eR50,xO(nmax),yO(nmax) DOUBLE PRECISION UPSF(nmax),pi,beta,value DOUBLE PRECISION alpha,xPSF(nmax),z(nmax),SB50,eSB50 DOUBLE PRECISION t(nmax),xp(nmax),yp(nmax),zp(nmax) DOUBLE PRECISION vx(nmax),vy(nmax),sigma(nmax),wk(nmax) DOUBLE PRECISION xsp(nmax),ysp(nmax),cfit(2) DOUBLE PRECISION xfc(nmax),yfc(nmax),sig(nmax),nsig(nmax) DOUBLE PRECISION convC(nmax),res(nmax) DOUBLE PRECISION OmSB50(nmax),OmR50(nmax),oR100,oI100 DOUBLE PRECISION FVEC(nmax),WA(nmax),DPMPAR DOUBLE PRECISION seta(4),sR50(4),schi(4) DOUBLE PRECISION chi(250),sSB50(4),sminim(4) DOUBLE PRECISION mini(nmax),yder(nmax),xder1,error INTEGER ider1 EXTERNAL FCN pi = 3.14159265358979323846 etamin = 0.3 etamax = 80 !10*etamax (eg: for etamax = 8 -> 80) opt_flag = 0 c = 0 scycle = 0 chi_flag = 7 calC = 25.D0 WRITE(*,*)' ' CALL getarg(1,filename) IF (filename.eq.' ') THEN WRITE(*,*)'Please insert name of file (x: arcsec; y: SB):' READ(*,*)filename ENDIF string = 'wc --l '//filename//' > tmp.dat' CALL system(string) OPEN(21,file='tmp.dat',status='old') READ(21,*)nlin0 CLOSE(21,status='delete') WRITE(*,*)'Filename: ',filename WRITE(*,*)'Number of rows: ',nlin0 string = 'wc --w '//filename//' > tmp.dat' CALL system(string) OPEN(21,file='tmp.dat',status='old') READ(21,*)ncol CLOSE(21,status='delete') ncol = ncol/nlin0 WRITE(*,*)'Number of columns: ',ncol CALL getarg(2,yPSF) IF (yPSF.eq.' ') THEN WRITE(*,*)'Want to take into account PFS conv. effects? (y/n)' READ(*,*)yPSF ENDIF WRITE(*,*)'PSF? ',yPSF IF (yPSF.eq.'y') THEN CALL getarg(3,MGU) IF (MGU.eq.' ') THEN WRITE(*,*)'Choose between Moffat (M) or Gaussian (G) or' WRITE(*,*)'to insert your own PSF (U) (x: arcsec; y: SB)' READ(*,*)MGU ENDIF IF (MGU.eq.'M') THEN WRITE(*,*)'Type of PSF: Moffat' ELSEIF (MGU.eq.'G') THEN WRITE(*,*)'Type of PSF: Gaussian' ELSEIF (MGU.eq.'U') THEN WRITE(*,*)'Type of PSF: User-defined' ENDIF IF (MGU.eq.'U') THEN CALL getarg(4,filenamePSF) IF (filenamePSF.eq.' ') THEN WRITE(*,*)'Please insert the name of the file' WRITE(*,*)'that contains the PSF (x: arcsec; y: SB):' READ(*,*)filenamePSF ENDIF string = 'wc --l '//filenamePSF//' > tmp.dat' CALL system(string) OPEN(21,file='tmp.dat',status='old') READ(21,*)nlinPSF CLOSE(21,status='delete') OPEN(90,file=filenamePSF,status='old') DO f = 1,nlinPSF READ(90,*)xPSF(f),UPSF(f) ENDDO CLOSE(90) DO f = 1,nlinPSF UPSF(f) = 10**((UPSF(f) - calC)/(-2.5D0)) ENDDO UPSF = UPSF / UPSF(1) WRITE(*,*)'File of PSF: ',filenamePSF WRITE(*,*)'Number of rows: ',nlinPSF ELSE CALL getarg(4,fw) IF (fw.eq.' ') THEN WRITE(*,*)'Insert FWHM (arcsec):' READ(*,*)FWHM ELSE READ(fw,*)FWHM ENDIF WRITE(*,'(A,F10.3)')' FWHM: ',FWHM ENDIF IF (MGU .eq. 'U') THEN cfit(1) = 4 cfit(2) = 4 lwa = nlinPSF*2 + 5*2 + nlinPSF val(2) = DSQRT(DPMPAR(1)) CALL LMDIF1(FCN,nlinPSF,2,cfit,FVEC,val(2),INFO,IWA, &WA,LWA,xPSF,UPSF) alpha = cfit(1) beta = cfit(2) FWHM = (2*alpha*sqrt(2**(1/beta)-1)) WRITE(*,'(A,F10.3)')' FWHM: ',FWHM WRITE(*,'(A,F10.3)')' beta: ',beta IF (alpha.le.0.or.beta.le.0) THEN WRITE(*,*)' ' WRITE(*,*)'WARNING:: THE PSF IS NOT PROPER!!' WRITE(*,*)' FITTING CANNOT PROCEED' STOP ELSEIF (alpha.ge.10.or.beta.ge.10) THEN WRITE(*,*)' ' WRITE(*,*)'WARNING:: THE PSF PARAMETERS ARE TOO HIGH!!' WRITE(*,*)' Fitting will proceed but its advised' WRITE(*,*)' to try different PSF' ENDIF ENDIF ENDIF IF (yPSF.eq.'y') THEN CALL getarg(5,fw) ELSE CALL getarg(3,fw) ENDIF IF (fw.eq.' ') THEN WRITE(*,*)'To compute distance dependent quantities insert' WRITE(*,*)'distance (Mpc) (not mandatory, you may give 0)' READ(*,*)distance IF (distance.ne.0) THEN WRITE(*,'(A,F10.3)')' Distance: ',distance ENDIF ELSE READ(fw,*)distance IF (distance.ne.0) THEN WRITE(*,'(A,F10.3)')' Distance: ',distance ENDIF ENDIF OPEN(99,file=filename,status='old') DO i = 1,nlin0 IF (ncol.eq.2) THEN READ(99,*)xi(i),yi(i) ELSEIF (ncol.eq.3) THEN READ(99,*)xi(i),yi(i),sig(i) ENDIF ENDDO npoints = nlin0 xstep(1) = (xi(10) - xi(3))/7 IF (npoints.gt.nmax) THEN WRITE(*,*)"WARNING:: TOO MANY POINTS" npoints = nmax ENDIF WRITE(*,*)" " WRITE(*,'(A,F10.3)')" Limiting radius: ",xi(npoints) WRITE(*,'(A,F10.3)')" Limiting SB: ",yi(npoints) rmin = xi(1) rmax = xi(npoints) !maximum obs. radius ! DETERMINE THE Nº OF EXTRAPOLATION POINTS jstep = npoints*4 IF (jstep.gt.2000) jstep=2000 IF (jstep.lt.500) jstep=500 val(1) = (rmax-rmin)/jstep DO i = 1,jstep xsp(i) = rmin + val(1)*dfloat(i-1) ENDDO xstep(2) = xsp(3) - xsp(2) !! SPLINE FOR SIGMA IF (maxval(sig).ne.0) THEN CALL ARCL2D(npoints,xi,sig,t,ier) j = 4*(npoints-1) CALL TSPSP(npoints,2,xi,sig,z,2,0,.FALSE.,.FALSE., &j,wk,t,xp,yp,zp,sigma,ier) CALL TSVAL1(npoints,xi,sig,yp,sigma,0,jstep,xsp, &nsig,IER) nsig = abs(nsig) ENDIF CALL MFlux(npoints,xi,yi,Rx,Ix,lum,calC,jstep,xsp,ysp) xO = xsp yO = ysp !! COMPUTE OBS DERIVATIVE a = 0 DO i = 2,jstep yder(i) = (yO(i+1) - yO(i))/(xO(i+1) - xO(i)) IF (i.gt.3.and.yder(i).gt.yder(i-1).and.a.eq.0) THEN a = 1 ider1 = i EXIT ENDIF ENDDO xder1 = xO(ider1) mag = -2.50D0*dlog10(lum) + calC !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! oR10 = Rx(1) oR20 = Rx(2) oR30 = Rx(3) oR40 = Rx(4) oR50 = Rx(5) oR60 = Rx(6) oR70 = Rx(7) oR80 = Rx(8) oR90 = Rx(9) oR100 = Rx(10) oI10 = Ix(1) oI20 = Ix(2) oI30 = Ix(3) oI40 = Ix(4) oI50 = Ix(5) oI60 = Ix(6) oI70 = Ix(7) oI80 = Ix(8) oI90 = Ix(9) oI100 = Ix(10) WRITE(*,*)' ' WRITE(*,*)'Observed Quantities:' WRITE(*,'(A,F10.3,A,F10.3)')' R40:',oR40,' :',oI40 WRITE(*,'(A,F10.3,A,F10.3)')' R50:',oR50,' :',oI50 WRITE(*,'(A,F10.3,A,F10.3)')' R60:',oR60,' :',oI60 WRITE(*,'(A,F10.3,A,F10.3)')' R70:',oR70,' :',oI70 WRITE(*,'(A,F10.3,A,F10.3)')' R80:',oR80,' :',oI80 WRITE(*,'(A,F10.3,A,F10.3)')' R90:',oR90,' :',oI90 WRITE(*,'(A,F10.3,A,F10.3)')' R100:',oR100,' :',oI100 WRITE(*,*)' ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 500 c = c + 1 g = jstep IF (yPSF.eq.'y') THEN k = minloc(abs(xO-3*FWHM),dim=1, mask = (xO .gt. 0.0000001)) IF (k.eq.g) THEN WRITE(*,*)' ' WRITE(*,*)'WARNING:: Cant assure PSF wont affect' WRITE(*,*)' collected data for Reff determination' WRITE(*,*)' ' k = minloc(abs(xO-2.5*FWHM),dim=1, mask = (xO .gt. 0.0000001)) IF (k.eq.g) k = minloc(abs(xO-FWHM),dim=1, &mask = (xO .gt. 0.0000001)) IF (k.eq.g) k = 1 ENDIF ELSE k = minloc(abs(xO-xstep(1)),dim=1, mask = (xO .gt. 0.0000001)) ENDIF !! REGRESSION TO FIX R50/SB50 FOR EACH ETA DO i = 1,etamax etaT(i) = etamin + (i-1)*0.1 CALL CompBeta(etaT(i),bn) xmod = xO**(1/etaT(i)) IF (ncol.eq.3) THEN IF (c.eq.1) THEN CALL LinearReg(k,g,xmod,yO,slo,b,eslo,eb) ELSE CALL LinearRegS(k,g,xmod,yO,nsig,slo,b,c,eslo,eb) ENDIF ELSE CALL LinearReg(k,g,xmod,yO,slo,b,eslo,eb) ENDIF OmSB50(i) = b + 2.5D0*bn/log(10.0D0) OmR50(i) = ((OmSB50(i)-b)/slo)**etaT(i) ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! run = 0 etalim(1) = 1 etalim(2) = 9 minim = 10000 50 n = 1 IF (yPSF.eq.'y') THEN run = run + 1 ! Prepare and convolve models DO g = etalim(1),etalim(2)!first run goes from etamin to etamax in steps of 1 i = g IF (minim(i).ne.10000) GOTO 10 IF (run.eq.1.and.g.gt.1) i = 10*g - 12 etaT(i) = etamin + (i-1)*0.1 CALL CompBeta(etaT(i),bn) R50 = OmR50(i) SB50 = OmSB50(i) IF (R50.eq.0.or.R50.gt.1000.or.R50.ne.R50) GOTO 10 DO f = 1,nmax/2 xmod(f) = (f-1)*xstep(1) ymod(f) = (2.5D0*bn/dlog(10.0D0)) ymod(f) = SB50 + ((xmod(f)/R50)** &(1.D0/etaT(i)) - 1.D0)*ymod(f) IF (xmod(f).ge.rmax+10) EXIT ENDDO lim(1) = f IF (lim(1).gt.(nmax-100)/2) lim(1) = (nmax-100)/2 lim(2) = lim(1)*2 - 1 DO f = 1,lim(2) IF (f.le.lim(1)) THEN xsym(f) = xmod(lim(1)-(f-1)) ysym(f) = ymod(lim(1)-(f-1)) ELSE xsym(f) = xmod((f-lim(1)+1)) ysym(f) = ymod((f-lim(1)+1)) ENDIF END DO DO f = 1,lim(2) ysymc(f) = 10**((ysym(f) - calC) / (-2.5D0)) ENDDO PSF = 0 ! CONSTRUCT PSF ACCORDING TO USER'S CHOICE IF (MGU .eq. 'M') THEN beta = 4.765 ! Trujillo et al. (2001) ELSEIF (MGU .eq. 'G') THEN beta = 1000000 ! Gaussian when beta -> inf ENDIF alpha = FWHM/(2.D0*sqrt(2.D0**(1.0D0/beta)-1.0D0)) DO f = 1,lim(2) PSF(f) = (1.0D0 + (xsym(f)/alpha)**2.D0)**(-1.D0*beta) ENDDO PSF = PSF/sum(PSF) ! MIDAS::: convC = 0 conv = 0 CALL CONVL1(ysymc,lim(2),PSF,lim(2),convC) DO f = lim(1),lim(2) conv(f-lim(1)+1) = -2.5D0*dlog10(convC(f)) + calC ENDDO !STOP THE CONV DATA DO f = 1,lim(1) IF (xmod(f).ge.rmax.or.conv(f).eq.0) THEN lim(1) = f-1 EXIT ELSE lim(1) = f ENDIF ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL MFlux(lim(1),xmod,conv,Rx,Ix,lum,calC,jstep,xsp,ysp) !! Xi2 not only for 1st point but for several initial points mini = 0 DO f = 1,1000 IF (xO(f).le.xder1) THEN val(1) = 10.0D0**((yO(f) - calC)/(-2.5D0)) val(2) = 10.0D0**((ysp(f) - calC)/(-2.5D0)) mini(f) = (val(1) - val(2))**2 ELSE EXIT ENDIF ENDDO error = sum(mini)/f minim(i) = minim(i) + error error = (oI40 - Ix(4))**2 minim(i) = minim(i) + error error = (oI50 - Ix(5))**2 minim(i) = minim(i) + error error = (oI60 - Ix(6))**2 minim(i) = minim(i) + error error = (oI70 - Ix(7))**2 minim(i) = minim(i) + error error = (oI80 - Ix(8))**2 minim(i) = minim(i) + error error = (oI90 - Ix(9))**2 minim(i) = minim(i) + error error = (oI100 - Ix(10))**2 minim(i) = minim(i) + error minim(i) = sqrt(minim(i)-10000)/8 + 1 ! GOODNESS OF FIT !!!!!!!!!!!!!!!!!!!! ! SPLINE CONVOLVED MODEL TO OBSERVATION CALL ARCL2D(lim(1),xmod,conv,t,ier) j = 4*(lim(1)-1) CALL TSPSP(lim(1),2,xmod,conv,z,2,0,.FALSE.,.FALSE., &j,wk,t,xp,yp,zp,sigma,ier) oy = 0 CALL TSVAL1(lim(1),xmod,conv,yp,sigma,0,jstep, &xO,oy,IER) ! COMPUTE DIFFERENCE BETWEEN MODEL AND OBSERVATION t = 0 DO f = 1,jstep t(f) = (yO(f) - oy(f))**2 ENDDO chi(i) = sqrt(sum(t)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (run.eq.1) n = 10 IF (i.gt.1.and.minim(i).gt.minim(i-n)) THEN m = m + 1 ELSE m = 0 ENDIF !! if 2 consecutive times minim increases no need to continue IF (m.eq.2.and.minim(i).gt.minim(i-n)) EXIT 10 ENDDO !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ELSE run = run + 1 DO g = etalim(1),etalim(2)!first run goes from etamin to etamax in steps of 1 i = g IF (minim(i).ne.10000) GOTO 11 IF (run.eq.1.and.g.gt.1) i = 10*g - 12 SB50 = OmSB50(i) R50 = OmR50(i) etaT(i) = etamin + (i-1)*0.1 IF (etaT(i).eq.1) etaT(i) = 1.000001 CALL CompBeta(etaT(i),bn) DO f = 1,nmax xmod(f) = (f-1)*xstep(1) ymod(f) = (2.5D0*bn/dlog(10.0D0)) ymod(f) = SB50 + ((xmod(f)/R50)** &(1.D0/etaT(i)) - 1.D0)*ymod(f) IF (xmod(f).ge.rmax) EXIT ENDDO IF (f.gt.(nmax)) f = nmax lim(1) = f ! MINIMIZATION FOR NON-CONV DATA:: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CALL MFlux(lim(1),xmod,ymod,Rx,Ix,lum,calC,jstep,xsp,ysp) error = (oI40 - Ix(4))**2 minim(i) = minim(i) + error error = (oI50 - Ix(5))**2 minim(i) = minim(i) + error error = (oI60 - Ix(6))**2 minim(i) = minim(i) + error error = (oI70 - Ix(7))**2 minim(i) = minim(i) + error error = (oI80 - Ix(8))**2 minim(i) = minim(i) + error error = (oI90 - Ix(9))**2 minim(i) = minim(i) + error error = (oI100 - Ix(10))**2 minim(i) = minim(i) + error minim(i) = sqrt(minim(i)-10000)/7 + 1 IF (run.eq.1) n = 10 IF (i.gt.1.and.minim(i).gt.minim(i-n)) THEN m = m + 1 ELSE m = 0 ENDIF ! COMPUTE DIFFERENCE BETWEEN MODEL AND OBSERVATION t = 0 DO f = 1,jstep t(f) = (yO(f) - ysp(f))**2 ENDDO chi(i) = sqrt(sum(t)) 11 ENDDO ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Find eta f = minloc(minim, dim=1, mask = (minim .gt. 0.00000000001)) !! SAVE ETA,R50,SB50 ESTIMATES TO ESTIMATE ERROR IF (run.eq.2) THEN seta(c) = etaT(f) sR50(c) = OmR50(f) sSB50(c) = OmSB50(f) sminim(c) = minim(f) schi(c) = chi(f) ENDIF IF (run.eq.1) THEN etalim(1) = f - 9 IF (etalim(1).lt.1) etalim(1) = 2 etalim(2) = f + 9 IF (etalim(2).gt.etamax) etalim(2) = etamax-1 GOTO 50 ENDIF IF (c.eq.1.and.ncol.eq.3) THEN WRITE(*,*)' :: Estimating uncertainties' scycle = 1 ENDIF IF (ncol.eq.2) THEN scycle = 0 g = 1 eta = seta(g) CALL CompBeta(eta,bn) R50 = sR50(g) SB50 = sSB50(g) chi(f) = schi(g) ENDIF IF (run.eq.2.and.scycle.eq.1) THEN IF (c.le.2) GOTO 500 IF (c.eq.3) GOTO 500 IF (c.eq.4) THEN g = minloc(sminim, dim=1) eta = seta(g) CALL CompBeta(eta,bn) R50 = sR50(g) SB50 = sSB50(g) chi(f) = schi(g) val(1) = sum(seta)/4 val(2) = sum(sR50)/4 val(3) = sum(sSB50)/4 val(4) = 0 val(5) = 0 val(6) = 0 DO i = 1,4 val(4) = val(4) + (seta(i) - val(1))**2 val(5) = val(5) + (sR50(i) - val(2))**2 val(6) = val(6) + (sSB50(i) - val(3))**2 ENDDO eeta = sqrt(val(4)/3) eR50 = sqrt(val(5)/3) eSB50 = sqrt(val(6)/3) ENDIF ENDIF IF (chi(f).gt.chi_flag) THEN IF (chi(f).gt.10) THEN opt_flag = 2 ELSEIF (chi(f).gt.chi_flag.and.chi(f).lt.10) THEN opt_flag = 1 ENDIF ENDIF DO i = 1,nmax xf(i) = (i - 1)*xstep(2) yf(i) = SB50 + (2.5*bn/log(10.0D0)) &*((xf(i)/R50)**(1.D0/eta)-1.D0) IF (yf(i).ge.30.or.i.gt.(nmax-100)/2) EXIT ENDDO lim(1) = i IF (yPSF.eq.'y') THEN DO i = 1,nmax xfc(i) = (i - 1)*xstep(1) yfc(i) = SB50 + (2.5*bn/log(10.0D0))*((xfc(i)/R50) &**(1.D0/eta)-1.D0) IF (yfc(i).ge.30.or.i.gt.(nmax-100)/2) EXIT ENDDO lim(3) = i C Create the convolved profile for output !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! lim(2) = lim(3)*2 - 1 DO i = 1,lim(2) IF (i.le.lim(3)) THEN xsym(i) = xfc(lim(3)-(i-1)) ysym(i) = yfc(lim(3)-(i-1)) ELSE xsym(i) = xfc((i-lim(3)+1)) ysym(i) = yfc((i-lim(3)+1)) ENDIF END DO DO i = 1,lim(2) ysymc(i) = 10**((ysym(i) - calC) / (-2.5D0)) ENDDO C Construct PSF according to user's choice: IF (MGU .eq. 'M') THEN beta = 4.765 ! Trujillo et al. (2001) ELSEIF (MGU .eq. 'G') THEN beta = 1000000 ! Gaussian when beta -> inf ENDIF alpha = FWHM/(2.D0*sqrt(2.D0**(1.0D0/beta)-1.0D0)) PSF = 0 DO i = 1,lim(2) PSF(i) = (1.D0 + (xsym(i)/alpha)**2.D0)**(-1.D0*beta) ENDDO PSF = PSF/sum(PSF) !! MIDAS :: convC = 0 conv = 0 CALL CONVL1(ysymc,lim(2),PSF,lim(2),convC) DO i = lim(3),lim(2) conv(i-lim(3)+1) = -2.50D0*dlog10(convC(i)) + calC ENDDO ENDIF !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IF (opt_flag.gt.1.and.yPSF.eq.'y') THEN WRITE(*,*)' ADVICE :: you may want to try a different PSF' WRITE(*,*)' ' ENDIF C Save model profile in ASCII file: f = index(filename,' ')-1 OPEN(82,file=filename(1:f)//'.SProf.dat',status='replace') DO i = 1,lim(1) IF (yPSF.eq.'y') THEN IF (i.le.lim(3)) THEN WRITE(82,*)xf(i),yf(i),xfc(i),conv(i) ELSE WRITE(82,*)xf(i),yf(i) ENDIF ELSE WRITE(82,*)xf(i),yf(i) ENDIF ENDDO CLOSE(82) WRITE(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' WRITE(*,*)' RESULT:' IF (ncol.eq.3) THEN WRITE(*,'(A,F10.1,A,F8.2)')' Sersic eta:',eta, &' +/-',eeta WRITE(*,'(A,F10.3,A,F8.2)')' Reff :',R50, &' +/-',eR50 ELSE WRITE(*,'(A,F10.1)')' Sersic eta:',eta WRITE(*,'(A,F10.3)')' Reff :',R50 ENDIF IF (distance.ne.0) THEN val(1) = 2.D0*pi*distance*1000.D0/(3600.D0*360.D0) val(2) = R50*val(1) val(3) = eR50*val(1) WRITE(*,'(A,F10.3,A,F8.2)')' Reff (kpc):', &val(2),' +/-',val(3) ENDIF IF (ncol.eq.3) THEN WRITE(*,'(A,F10.3,A,F8.2)')' SBReff :', &SB50,' +/-',eSB50 ELSE WRITE(*,'(A,F10.3)')' SBReff :',SB50 ENDIF WRITE(*,'(A,F10.3)')' AppM :',mag IF (distance.ne.0) THEN magAbs = mag - (calC + (5.D0*log10(distance))) + 5 WRITE(*,'(A,F10.3)')' AbsM :',magAbs ENDIF WRITE(*,'(A,F10.3)')' SB(R=0) :',yf(1) WRITE(*,*)' ' WRITE(*,'(A,A,A)')' Modeled profile in file ',filename(1:f)// &'.Ser_prof.dat' WRITE(*,'(A,A,A)')' Fit results in file ',filename(1:f)// &'.Results.dat' WRITE(*,*)' ' WRITE(*,*)'Modeled Quantities:' WRITE(*,'(A,F10.3,A,F10.3)')' R40:',Rx(4),' :',Ix(4) WRITE(*,'(A,F10.3,A,F10.3)')' R50:',Rx(5),' :',Ix(5) WRITE(*,'(A,F10.3,A,F10.3)')' R60:',Rx(6),' :',Ix(6) WRITE(*,'(A,F10.3,A,F10.3)')' R70:',Rx(7),' :',Ix(7) WRITE(*,'(A,F10.3,A,F10.3)')' R80:',Rx(8),' :',Ix(8) WRITE(*,'(A,F10.3,A,F10.3)')' R90:',Rx(9),' :',Ix(9) WRITE(*,'(A,F10.3,A,F10.3)')' R100:',Rx(10), &' :',Ix(10) WRITE(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' OPEN(99,file=filename(1:f)//'.Results.dat',status='replace') IF (distance.eq.0) THEN WRITE(99,'(F10.1,F10.3,F10.3,F10.3,F10.3)')eta,R50,SB50, &mag,yf(1) ELSE WRITE(99,'(F10.1,F10.3,F10.3,F10.3,F10.3,F10.3,F10.3)')eta, &R50,val(2),SB50,mag,magAbs,yf(1) ENDIF CLOSE(99) C Plot the result with PGPLOT ------------------------------------- CALL PGBEGIN(0,filename(1:f)//'.ps/cps',1,1) C -------- Define new colors ----------------- CALL PGSCR(22,0.97,0.97,0.97) ! open gray CALL PGSCR(32,0.90,0.90,0.90) ! open gray 85% CALL PGSCR(33,0.789,0.0436,0.0833) ! OpenOffice red/chart CALL PGSCR(34,0.0436,0.54365,0.8333) ! OpenOffice blue/chart CALL PGSCR(35,0.7411,0.776,0.1098) ! my OpenGreen 1 CALL PGSCR(36,0.579,0.769,0.675) ! dark green CALL PGSCR(30,0.65,0.65,0.65) ! mid gray CALL PGSCR(23,0.2,0.3,0.4) ! dark gray CALL PGSCR(18,0.88,0.94,0.94) ! Himmelsblau CALL PGSCR(21,0.246,0.62,1.0) ! open blue (11) CALL PGSCR(19,0.349,0.539,0.678) ! blue jeans CALL PGSCR(25,0.091,0.1548,0.254) ! metallic blue CALL PGSCR(20,0.2941,1.0,0.8) ! open green CALL PGSCR(17,0.62,0.73,0.68) ! chaki CALL PGSCR(16,0.937,0.4941,0.0) ! dark orange CALL PGSCR(24,0.4,0.095,0.119) ! metallic red CALL PGSCR(26,0.73,0.73,0.73) ! Mittelgrau CALL PGSCR(27,0.15,0.3,0.15) ! very dark grey/green CALL PGSCR(50,0.992,0.5216,0.5216) ! salmon CALL PGSCR(51,0.6196,0.8078,1.0) ! OO hellblau CALL PGSCR(52,0.50,0.00,1.00) ! purple C -------- Background color ----------------- CALL PGSCR(0,1.0,1.0,1.0) ! white background CALL PGSCR(1,0.0,0.0,0.0) ! define black C --------------------------------------------------- CALL PGSCF(2) ! ifont CALL PGSLW(4) ! line width I*2 CALL PGSCH(0.67) ! ssize C Legend: CALL PGSCI(0) CALL PGSLW(1) CALL PGSVP(0.65,0.86,0.8,0.88) CALL PGSWIN(0.0, 1.0, 0.0, 1.0) CALL PGBOX('BCNTS',0.0,0,'BCNTSV',0.0,0) CALL PGSCI(1) CALL PGSLW(2) CALL PGPT(1, 0.1, 0.9, -4) CALL PGTEXT(0.2, 0.86, '\fnInput Data Points') vx(1) = 0.05 vx(2) = 0.15 vy(1) = 0.65 vy(2) = 0.65 CALL PGSLW(5) CALL PGSCI(33) CALL PGLINE(2, real(vx), real(vy)) CALL PGSCI(1) CALL PGSLW(2) CALL PGTEXT(0.2, 0.61, '\fnRetrieved Sersic Profile') IF (yPSF.eq.'y') THEN CALL PGSCI(19) CALL PGSLW(5) vy(1) = 0.40 vy(2) = 0.40 CALL PGLINE(2, real(vx), real(vy)) CALL PGSCI(1) CALL PGSLW(2) CALL PGTEXT(0.2, 0.36, '\fnConvolved Sersic Profile') ENDIF C --------------------------------------- C Limits of main plot ymax = yf(1) - 1.5 value = xi(nlin0)*1.2 i = minloc(abs(xf-value),dim=1, mask = (xf .gt. 0.0000001)) xmax = xf(i) ymin = yf(i) CALL PGSLW(2) ! line thickness CALL PGSVP(0.1,0.9,0.23,0.93) CALL PGSWIN(0., xmax, ymin, ymax) CALL PGSCI(1) CALL PGBOX('BCTS',0.0,0,'BCNTSV',0.0,0) f = index(filename,' ')-1 IF (yPSF.eq.'y') THEN WRITE(fw,'(F4.2)')FWHM WRITE(string,'(F4.2)')beta IF (MGU.eq.'U') THEN k = index(filenamePSF,' ')-1 CALL PGLAB('', '', &filename(1:f)//' PSF: '//filenamePSF(1:k)// &' FWHM :'//fw//' \gb : '//string) ELSEIF (MGU.eq.'M') THEN CALL PGLAB('', '', &filename(1:f)//' FWHM :'//fw// &' \gb : '//string) ELSEIF (MGU.eq.'G') THEN CALL PGLAB('', '', &filename(1:f)//' FWHM :'//fw// &' \gb : \(0766)') ENDIF ELSE CALL PGLAB('', '', &filename(1:f)) ENDIF CALL PGMTXT ('L', 2.4, 0.4, 1, '\gm (mag/arcsec\u2)') CALL PGPT(nlin0, real(xi), real(yi), -4) IF (ncol.eq.3) THEN WRITE(string,'(F6.1,A,F6.2)')eta,' +/- ',eeta ELSE WRITE(string,'(F6.1)')eta ENDIF CALL PGMTXT ('T', -1.5, 0.05, 1, ' \gy : '//string) IF (ncol.eq.3) THEN WRITE(string,'(F6.2,A,F6.2)')real(R50),' +/- ',eR50 ELSE WRITE(string,'(F6.2)')real(R50) ENDIF CALL PGMTXT ('T', -2.7, 0.05, 1, 'Reff : '//string) IF (ncol.eq.3) THEN WRITE(string,'(F6.2,A,F6.2)')real(SB50),' +/- ',eSB50 ELSE WRITE(string,'(F6.2)')real(SB50) ENDIF CALL PGMTXT ('T', -3.9, 0.05, 1, '\gmeff : '//string) C Watermark CALL PGSLW(1) CALL PGSCI(32) CALL PGMTXT ('B', 10.0, 1, 1, '\(2380)ris') CALL PGSLW(4) C Method to estimate R50 & mu50 -- color of the inner circle C - blue : simple linear regression C - purple : w = 1/sqrt(sigma) C - orange : w = 1/sigma C - red : w = 1/sigma² IF (g.eq.1) THEN CALL PGSCI(19) !blue CALL PGMTXT ('T', 1.85, 0.97, 1, '\(0828)') ELSEIF (g.eq.2) THEN CALL PGSCI(52) !purple CALL PGMTXT ('T', 1.85, 0.97, 1, '\(0828)') ELSEIF (g.eq.3) THEN CALL PGSCI(16) !orange CALL PGMTXT ('T', 1.85, 0.97, 1, '\(0828)') ELSEIF (g.eq.4) THEN CALL PGSCI(33) !red CALL PGMTXT ('T', 1.85, 0.97, 1, '\(0828)') ENDIF C Goodness of the fit - color of the outer circle C - blue : |mean(Chi) + st(Chi)|/2 < 1 - good C - purple : 1 < |mean(Chi) + st(Chi)|/2 < 2 - moderate C - red : |mean(Chi) + st(Chi)|/2 > 2 - bad IF (opt_flag.eq.0) THEN CALL PGSCI(19) CALL PGMTXT ('T', 1.85, 0.963, 1, '\(0905)') ELSEIF (opt_flag.eq.2) THEN CALL PGSCI(33) CALL PGMTXT ('T', 1.85, 0.963, 1, '\(0905)') ELSEIF (opt_flag.eq.1) THEN CALL PGSCI(52) CALL PGMTXT ('T', 1.85, 0.963, 1, '\(0905)') ENDIF C -------------------------------------------------------------- CALL PGSLS(1) !-> solid line CALL PGSCI(1) CALL PGSLW(2) C Observational points: IF (ncol.eq.3) THEN CALL PGERRY(nlin0, real(xi), real(yi-sig), real(yi+sig), 1.) ENDIF C Model: CALL PGSCI(33) CALL PGSLW(5) CALL PGLINE(lim(1), real(xf), real(yf)) C Convolved Model: IF (yPSF.eq.'y') THEN CALL PGSCI(19) CALL PGSLW(5) CALL PGLINE(lim(3), real(xfc), real(conv)) ENDIF C Subplot: CALL PGSLW(2) CALL PGSVP(0.45,0.88,0.5,0.9) IF (yPSF.eq.'y'.and.MGU.eq.'U') THEN ymax = UPSF(1) + (UPSF(1))*0.1 i = nlinPSF CALL PGSWIN(0.,real(xPSF(i)),real(UPSF(i)),ymax) CALL PGSCI(1) CALL PGBOX('BCNTS',0.0,0,'BCNTSV',0.0,0) CALL PGMTXT ('T', 0.3, 0., 1, '\fnUser Defined PSF:') C Observed PSF: CALL PGSCI(1) CALL PGPT(nlinPSF, real(xPSF), real(UPSF), -4) C Model PSF: CALL PGSCI(33) CALL Moffat(nlinPSF,xPSF,cfit,Mof) CALL PGLINE(nlinPSF, real(xPSF), real(Mof)) ELSEIF (eta.ge.1.5) THEN C Zoom in the 1st 10% of input data: ymax = yf(1) - 1 i = minloc(abs(xi-0.1*(xi(nlin0))),dim=1, mask = &(xi .gt. 0.0000001)) CALL PGSWIN(0.,real(xi(i)),real(yi(i)),ymax) CALL PGSCI(1) CALL PGBOX('BCNTS',0.0,0,'BCNTSV',0.0,0) CALL PGMTXT ('T', 0.3, 0., 1, '\fnZoom at Inner Profile:') CALL PGSCI(1) CALL PGSLW(2) CALL PGPT(nlin0, real(xi), real(yi), -4) CALL PGSLW(5) CALL PGSCI(33) CALL PGLINE(lim(1), real(xf), real(yf)) IF (yPSF.eq.'y') THEN CALL PGSCI(19) CALL PGLINE(lim(3), real(xfc), real(conv)) ENDIF ENDIF C Residuals: IF (yPSF.eq.'y') THEN CALL ARCL2D(lim(3),xfc,conv,t,ier) j = 4*(lim(3)-1) CALL TSPSP(lim(3),2,xfc,conv,z,2,0,.FALSE.,.FALSE., &j,wk,t,xp,yp,zp,sigma,ier) oy = 0 CALL TSVAL1(lim(3),xfc,conv,yp,sigma,0,jstep, &xO,oy,IER) ELSE CALL ARCL2D(lim(1),xf,yf,t,ier) j = 4*(lim(1)-1) CALL TSPSP(lim(1),2,xf,yf,z,2,0,.FALSE.,.FALSE., &j,wk,t,xp,yp,zp,sigma,ier) oy = 0 CALL TSVAL1(lim(1),xf,yf,yp,sigma,0,jstep, &xO,oy,IER) ENDIF res = yO - oy CALL PGSLW(2) CALL PGSVP(0.1,0.9,0.12,0.22) CALL PGSWIN(0.0, xmax, -1.0, 1.0) CALL PGSCI(1) CALL PGBOX('BCNTS',0.0,0,'BCNTSV',0.0,0) CALL PGLAB('Radius (arcsec)', '', '') CALL PGSLW(5) CALL PGLINE(jstep, real(xsp), real(res)) CALL PGSLW(2) CALL PGMTXT ('L', 3.1, -0.2, 1, 'Residuals (mag)') C 0 line: res = 0 CALL PGSLS(5) CALL PGSLW(2) CALL PGLINE(lim(1), real(xf), real(res)) CALL PGEND END C ====================================================================================== SUBROUTINE MFlux(n,x,y,Rx,Ix,lum,calC,np,xsp,ysp) INTEGER i,f,n,ier PARAMETER (nmax=25000) REAL calC DOUBLE PRECISION x(n),y(n),pi,yc(n),wk(nmax) DOUBLE PRECISION t(nmax),xp(nmax),yp(nmax),zp(nmax),z(n) DOUBLE PRECISION xsp(nmax),ysp(nmax),sig(nmax),lum DOUBLE PRECISION tsintl,ratio(nmax),intens(nmax) DOUBLE PRECISION Rx(10),Ix(10) pi = 3.14159265358979323846 DO i = 1,n yc(i) = 2.0D0*pi*x(i)*(10.0D0**((y(i) - calC)/(-2.5D0))) END DO CALL ARCL2D(n,x,yc,t,ier) f = 4*(n-1) CALL TSPSP(n,2,x,yc,z,2,0,.FALSE.,.FALSE., &f,wk,t,xp,yp,zp,sig,ier) CALL TSVAL1(n,x,yc,yp,sig,0,np,xsp,ysp,IER) lum = TSINTL(xsp(1),xsp(np),np,xsp,ysp,yp,sig,ier) DO i = 1,np intens(i) = TSINTL(xsp(1),xsp(i),np,xsp,ysp,yp,sig,ier) ENDDO ratio = intens/lum yp = ysp f = minloc(abs(ratio - 0.1), dim=1, mask = (ratio .gt. 0.0000001)) Rx(1) = xsp(f) f = minloc(abs(ratio - 0.2), dim=1, mask = (ratio .gt. 0.0000001)) Rx(2) = xsp(f) f = minloc(abs(ratio - 0.3), dim=1, mask = (ratio .gt. 0.0000001)) Rx(3) = xsp(f) f = minloc(abs(ratio - 0.4), dim=1, mask = (ratio .gt. 0.0000001)) Rx(4) = xsp(f) f = minloc(abs(ratio - 0.5), dim=1, mask = (ratio .gt. 0.0000001)) Rx(5) = xsp(f) f = minloc(abs(ratio - 0.6), dim=1, mask = (ratio .gt. 0.0000001)) Rx(6) = xsp(f) f = minloc(abs(ratio - 0.7), dim=1, mask = (ratio .gt. 0.0000001)) Rx(7) = xsp(f) f = minloc(abs(ratio - 0.8), dim=1, mask = (ratio .gt. 0.0000001)) Rx(8) = xsp(f) f = minloc(abs(ratio - 0.9), dim=1, mask = (ratio .gt. 0.0000001)) Rx(9) = xsp(f) f = minloc(abs(ratio - 1), dim=1, mask = (ratio .gt. 0.0000001)) Rx(10) = xsp(f) ! MEAN INTENSITY: Ix(1) = 0.1*lum/(pi*Rx(1)**2) Ix(2) = 0.2*lum/(pi*Rx(2)**2) Ix(3) = 0.3*lum/(pi*Rx(3)**2) Ix(4) = 0.4*lum/(pi*Rx(4)**2) Ix(5) = 0.5*lum/(pi*Rx(5)**2) Ix(6) = 0.6*lum/(pi*Rx(6)**2) Ix(7) = 0.7*lum/(pi*Rx(7)**2) Ix(8) = 0.8*lum/(pi*Rx(8)**2) Ix(9) = 0.9*lum/(pi*Rx(9)**2) Ix(10) = lum/(pi*Rx(10)**2) ysp(1) = y(1) DO i = 2,np ysp(i) = (-2.5D0)*dlog10(yp(i)/(2.0D0*pi*xsp(i))) + calC ENDDO END C ====================================================================================== SUBROUTINE CONVL1(F,NF,PSF,NP,C) C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C subroutine CONVL1 version 1.00 880916 C P. GROSBOL ESO - GARCHING C K. Banse (streamlined for MIDAS) C C.KEYWORDS: C convolution C C.PURPOSE: C convoles two one dimensional arrays with each other C C.ALGORITHM: C a one dimensional array is convoled with a one dimensional C point spread function. outside the limits of the function C array the first/last value is taken. C C.INPUT/OUTPUT: C the subroutine is called as CONVL1(F,NF,PSF,NP,C) C C F(NF) : function array (INPUT) C NF : dimension of function array (INPUT) C PSF(NP) : point spread function array (INPUT) C NP : dimension of point spread function (INPUT) C C(NF) : convolved result array (OUTPUT) C C NOTE : the result is added to the output array 'C'. C C--------------------------------------------------------------- C IMPLICIT NONE C INTEGER NF,NP INTEGER I,K,NL1,NL2,NL3 INTEGER NOS,NR C DOUBLE PRECISION F(NF),C(NF),PSF(NP) DOUBLE PRECISION FACT,R C C initialize NOS = NP/2 NL1 = NOS + 1 NL2 = NF - NOS NL3 = NL2 + 1 C C if NP = 1, we only use middle section IF (NOS.EQ.0) THEN DO 1000, I=NL1,NL2 NR = I - NOS R = C(I) DO 400, K=1,NP R = R + PSF(K)*F(NR) NR = NR + 1 400 CONTINUE C(I) = R 1000 CONTINUE RETURN ENDIF C C for NP > 1, handle 1. section, where values before F(1) would be needed FACT = F(1) DO 2000, I=1,NOS R = C(I) DO 1400 K=1,1+NOS-I R = R + PSF(K)*FACT 1400 CONTINUE NR = 1 C DO 1600, K=2+NOS-I,NP R = R + PSF(K)*F(NR) NR = NR + 1 1600 CONTINUE C(I) = R 2000 CONTINUE C C in the middle no problems DO 3000 I=NL1,NL2 NR = I - NOS R = C(I) DO 2400, K=1,NP R = R + PSF(K)*F(NR) NR = NR + 1 2400 CONTINUE C(I) = R 3000 CONTINUE C C last section, values after F(NF) would be needed FACT = F(NF) DO 4000, I=NL3,NF NR = I - NOS R = C(I) DO 3300, K=1,NF+NOS-I R = R + PSF(K)*F(NR) NR = NR + 1 3300 CONTINUE C DO 3600, K=NF+NOS-I+1,NP R = R + PSF(K)*FACT 3600 CONTINUE C(I) = R 4000 CONTINUE C RETURN END C ====================================================================================== SUBROUTINE Moffat(n,x,c,Mof) IMPLICIT NONE INTEGER i,n DOUBLE PRECISION Mof(n),x(n),c(2),pi pi = 3.14159265358979323846 DO i = 1,n Mof(i) = (1.0D0 + (x(i)/c(1))**2.D0)**(-1.D0*c(2)) ENDDO END C ====================================================================================== ! Ciotti & Bertin (1999): ! http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1999A%26A...352..447C&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf SUBROUTINE CompBeta(eta,bn) REAL eta DOUBLE PRECISION bn bn = 2*eta - 1/3 + 4/(405*eta) + 46/(25515*eta**2) + &131/(1148175*eta**3) - 7.1510122958919723e-05/eta**4 END C -------------------------------------------------------------------------------------------------------------- SUBROUTINE LinearReg(lim1,lim2,ix,iy,slo,b,eslo,eb) !https://en.wikipedia.org/wiki/Simple_linear_regression INTEGER i,n,lim1,lim2 PARAMETER (nmax=25000) DOUBLE PRECISION ix(nmax),iy(nmax) DOUBLE PRECISION sum1(nmax),sum2(nmax),meanx,meany DOUBLE PRECISION sum3,sum4,sum5,sum6 DOUBLE PRECISION slo,eslo,b,eb,se n = lim2 - lim1 + 1 meanx = sum(ix(lim1:lim2))/n meany = sum(iy(lim1:lim2))/n DO i = lim1,lim2 sum1(i) = (ix(i) - meanx)*(iy(i) - meany) sum2(i) = (ix(i) - meanx)**2 ENDDO sum3 = sum(ix(lim1:lim2)) ! S(x) sum4 = sum(iy(lim1:lim2)) ! S(y) sum5 = sum(ix(lim1:lim2)**2) ! S(x²) sum6 = sum(iy(lim1:lim2)**2) ! S(y²) slo = sum(sum1(lim1:lim2))/sum(sum2(lim1:lim2)) b = meany - slo*meanx ! errors: se = n*(n-2) se = 1/se se = se * ( n*sum6 - sum4**2 - slo**2 * &(n*sum5 - sum3**2)) eslo = sqrt( n*se**2 / (n*sum5 - sum3**2) ) eb = sqrt( eslo**2 * sum5 / n ) END C -------------------------------------------------------------------------------------------------------------- SUBROUTINE LinearRegS(lim1,lim2,ix,iy,ysig,slo,b,h,eslo,eb) ! Linear Regression taking y errors into account !(https://www.originlab.com/doc/Origin-Help/FIt-with-Err-Weight) ! formula: ! (knoble_linfit.ps -> www.oc.nps.edu/~bird/oc3030_online/project/knoble_linfit.ps) ! ! http://www.lithoguru.com/scientist/statistics/Lecture28.pdf (WLR.pdf): ! w = 1/sig**2 INTEGER i,h,lim1,lim2,n PARAMETER (nmax=25000) DOUBLE PRECISION ix(nmax),iy(nmax),ysig(nmax) DOUBLE PRECISION w(nmax),sum1(nmax),sum2(nmax),sum3(nmax) DOUBLE PRECISION slo,eslo,b,eb,sum4(nmax),error(nmax) DOUBLE PRECISION val(10),sum5(nmax),sum6(nmax),se,xm DO i = lim1,lim2 IF (ysig(i).eq.0) ysig(i) = 1E-20 IF (h.eq.2) THEN w(i) = 1/ysig(i) ELSEIF (h.eq.3) THEN w(i) = 1/sqrt(ysig(i)) ELSEIF (h.eq.4) THEN w(i) = 1/ysig(i)**2 ENDIF sum1(i) = w(i)*ix(i) sum2(i) = w(i)*iy(i) sum3(i) = w(i)*ix(i)**2 sum4(i) = w(i)*iy(i)**2 sum5(i) = w(i)*ix(i)*iy(i) ENDDO val(1) = sum(w(lim1:lim2)) ! S(w) val(2) = sum(sum1(lim1:lim2)) ! S(w*x) val(3) = sum(sum2(lim1:lim2)) ! S(w*y) val(4) = sum(sum3(lim1:lim2)) ! S(w*x^2) val(5) = sum(sum4(lim1:lim2)) ! S(w*y^2) val(6) = sum(sum5(lim1:lim2)) ! S(w*x*y) val(7) = val(4)*val(1) - val(2)**2 ! S(w*x^2) * S(w) - S(w*x)^2 slo = (1/val(7))*(val(6)*val(1) - val(2)*val(3)) b = (1/val(7))*(val(4)*val(3) - val(6)*val(2)) ! errors: n = lim2 - lim1 + 1 xm = val(2)/val(1) DO i = lim1,lim2 error(i) = (sqrt(w(i))*(iy(i) - (b + slo*ix(i))))**2 sum6(i) = w(i)*(ix(i) - xm)**2 ENDDO val(8) = sum(error(lim1:lim2)) se = sqrt(val(8)/(n - 2)) val(9) = sum(sum6(lim1:lim2)) eslo = se/sqrt(val(9)) eb = sqrt(se**2/val(1) + (eslo*xm)**2) END C ---------------------------------------------------------- SUBROUTINE TSPSP (N,ND,X,Y,Z,NCD,IENDC,PER,UNIFRM, . LWK, WK, T,XP,YP,ZP,SIGMA,IER) INTEGER N, ND, NCD, IENDC, LWK, IER LOGICAL PER, UNIFRM DOUBLE PRECISION X(N), Y(N), Z(N), WK(LWK), T(N), . XP(N), YP(N), ZP(N), SIGMA(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 07/08/92 C C This subroutine computes a set of values which define a C parametric planar curve C(t) = (H1(t),H2(t)) or space C curve C(t) = (H1(t),H2(t),H3(t)) whose components are Her- C mite interpolatory tension splines. The output values C consist of parameters (knots) T computed by ARCL2D or C ARCL3D, knot derivative values XP, YP, (and ZP) computed C by Subroutine YPC1, YPC1P, YPC2, or YPC2P, and tension C factors SIGMA computed by Subroutine SIGS (unless UNIFRM = C TRUE). C C Refer to Subroutine TSPSP for an alternative method of C computing tension factors in the case of a planar curve. C C The tension splines may be evaluated by Subroutine C TSVAL2 (or TSVAL3) or Functions HVAL (values), HPVAL C (first derivatives), HPPVAL (second derivatives), and C TSINTL (integrals). C C On input: C C N = Number of knots and data points. N .GE. 2 and C N .GE. 3 if PER = TRUE. C C ND = Number of dimensions: C ND = 2 if a planar curve is to be constructed. C ND = 3 if a space curve is to be constructed. C C X,Y,Z = Arrays of length N containing the Cartesian C coordinates of an ordered sequence of data C points C(I), I = 1 to N, such that C(I) .NE. C C(I+1). C(t) is constrained to pass through C these points. Z is an unused dummy parame- C ter if ND = 2. In the case of a closed curve C (PER = TRUE), the first and last points C should coincide. In this case, X(N), Y(N), C (and Z(N)) are set to X(1), Y(1), (and Z(1)) C if NCD = 1, but not altered if NCD = 2. C C NCD = Number of continuous derivatives at the knots. C NCD = 1 or NCD = 2. If NCD = 1, XP, YP, (and C ZP) are computed by local monotonicity- C constrained quadratic fits. Otherwise, a C linear system is solved for the derivative C values which result in second derivative con- C tinuity. Unless UNIFRM = FALSE, this requires C iterating on calls to YPC2 or YPC2P and calls C to SIGS, and generally results in more nonzero C tension factors (hence more expensive evalua- C tion). C C IENDC = End condition indicator for NCD = 2 and PER C = FALSE (or dummy parameter otherwise): C IENDC = 0 if XP(1), XP(N), YP(1), YP(N) (and C ZP(1) and ZP(N)) are to be com- C puted by monotonicity-constrained C parabolic fits (YPC1). C IENDC = 1 if the first derivatives of H1 at C the left and right endpoints are C user-specified in XP(1) and XP(N), C respectively, the first deriva- C tives of H2 at the ends are C specified in YP(1) and YP(N), and, C if ND = 3, the first derivatives C of H3 are specified in ZP(1) and C ZP(N). C IENDC = 2 if the second derivatives of H1, C H2, (and H3) at the endpoints are C user-specified in XP(1), XP(N), C YP(1), YP(N), (ZP(1), and ZP(N)). C IENDC = 3 if the end conditions are to be C computed by Subroutine ENDSLP and C vary with SIGMA(1) and SIGMA(N-1). C C PER = Logical variable with value TRUE if and only C a closed curve is to be constructed -- H1(t), C H2(t), (and H3(t)) are to be periodic func- C tions with period T(N)-T(1), where T(1) and C T(N) are the parameter values associated with C the first and last data points. It is assumed C in this case that X(N) = X(1), Y(N) = Y(1) C and, if ND = 3, Z(N) = Z(1), and, on output, C XP(N) = XP(1), YP(N) = YP(1), (and ZP(N) = C ZP(1) if ND = 3). C C UNIFRM = Logical variable with value TRUE if and C only if constant (uniform) tension is to be C used. The tension factor must be input in C SIGMA(1) in this case and must be in the C range 0 to 85. If SIGMA(1) = 0, H(t) is C piecewise cubic (a cubic spline if NCD = C 2), and as SIGMA increases, H(t) approaches C the piecewise linear interpolant, where H C is H1, H2, or H3. If UNIFRM = FALSE, C tension factors are chosen (by SIGS) to C preserve local monotonicity and convexity C of the data. This often improves the C appearance of the curve over the piecewise C cubic fitting functions. C C LWK = Length of work space WK: no work space is C needed if NCD = 1; at least N-1 locations C are required if NCD = 2; another N-1 locations C are required if PER = TRUE; and an additional C ND*(N-1) locations are required for the con- C vergence test if SIGS is called (UNIFRM = C FALSE): C If NCD=1 then LWK = 0 (not tested). C If NCD=2 then C C LWK GE N-1 if PER=FALSE, UNIFRM=TRUE C LWK GE 2N-2 if PER=TRUE, UNIFRM=TRUE C LWK GE (ND+1)*(N-1) if PER=FALSE, UNIFRM=FALSE C LWK GE (ND+2)*(N-1) if PER=TRUE, UNIFRM=FALSE C C The above parameters, except possibly X(N), Y(N), and C Z(N), are not altered by this routine. C C WK = Array of length .GE. LWK to be used as tempor- C ary work space. C C T = Array of length .GE. N. C C XP = Array of length .GE. N containing end condition C values in positions 1 and N if NCD = 2 and C IENDC = 1 or IENDC = 2. C C YP = Array of length .GE. N containing end condition C values in positions 1 and N if NCD = 2 and C IENDC = 1 or IENDC = 2. C C ZP = Dummy argument if ND=2, or, if ND=3, array of C length .GE. N containing end condition values C in positions 1 and N if NCD = 2 and IENDC = 1 C or IENDC = 2. C C SIGMA = Array of length .GE. N-1 containing a ten- C sion factor (0 to 85) in the first position C if UNIFRM = TRUE. C C On output: C C WK = Array containing convergence parameters in the C first two locations if IER > 0 (NCD = 2 and C UNIFRM = FALSE): C WK(1) = Maximum relative change in a component C of XP, YP, or ZP on the last iteration. C WK(2) = Maximum relative change in a component C of SIGMA on the last iteration. C C T = Array containing parameter values computed by C Subroutine ARCL2D or ARCL3D unless -4 < IER < 0. C T is only partially defined if IER = -4. C C XP = Array containing derivatives of H1 at the C knots. XP is not altered if -5 < IER < 0, C and XP is only partially defined if IER = -6. C C YP = Array containing derivatives of H2 at the C knots. YP is not altered if -5 < IER < 0, C and YP is only partially defined if IER = -6. C C ZP = Array containing derivatives of H3 at the knots C if ND=3. ZP is not altered if -5 < IER < 0, C and ZP is only partially defined if IER = -6. C C SIGMA = Array containing tension factors. SIGMA(I) C is associated with interval (T(I),T(I+1)) C for I = 1,...,N-1. SIGMA is not altered if C -5 < IER < 0 (unless IENDC is invalid), and C SIGMA is constant (not optimal) if IER = -6 C or IENDC (if used) is invalid. C C IER = Error indicator or iteration count: C IER = IC .GE. 0 if no errors were encountered C and IC calls to SIGS and IC+1 calls C to YPC1, YPC1P, YPC2 or YPC2P were C employed. (IC = 0 if NCD = 1). C IER = -1 if N, NCD, or IENDC is outside its C valid range. C IER = -2 if LWK is too small. C IER = -3 if UNIFRM = TRUE and SIGMA(1) is out- C side its valid range. C IER = -4 if a pair of adjacent data points C coincide: X(I) = X(I+1), Y(I) = C Y(I+1), (and Z(I) = Z(I+1)) for some C I in the range 1 to N-1. C IER = -6 if invalid knots T were returned by C ARCL2D or ARCL3D. This should not C occur. C C Modules required by TSPSP: ARCL2D, ARCL3D, ENDSLP, SIGS, C SNHCSH, STORE, YPCOEF, YPC1, C YPC1P, YPC2, YPC2P C C Intrinsic functions called by TSPSP: ABS, MAX C C*********************************************************** C INTEGER I, ICNT, IERR, ITER, IW1, MAXIT, N2M2, NM1, NN LOGICAL SCURV DOUBLE PRECISION DSMAX, DYP, DYPTOL, EX, EY, EZ, SBIG, . SIG, STOL, XP1, XPN, YP1, YPN, ZP1, . ZPN C DATA SBIG/85.D0/ C C Convergence parameters: C C STOL = Absolute tolerance for SIGS C MAXIT = Maximum number of YPC2/SIGS iterations C DYPTOL = Bound on the maximum relative change in a com- C ponent of XP, YP, or ZP defining convergence C of the YPC2/SIGS iteration when NCD = 2 and C UNIFRM = FALSE C DATA STOL/0.D0/, MAXIT/99/, DYPTOL/.01D0/ C C Test for invalid input parameters N, NCD, or LWK. C NN = N NM1 = NN - 1 N2M2 = 2*NM1 IF (NN .LT. 2 .OR. (PER .AND. NN .LT. 3) .OR. . NCD .LT. 1 .OR. NCD .GT. 2) GO TO 11 IF (UNIFRM) THEN IF ( NCD .EQ. 2 .AND. (LWK .LT. NM1 .OR. . (PER .AND. LWK .LT. N2M2)) ) GO TO 12 SIG = SIGMA(1) IF (SIG .LT. 0.D0 .OR. SIG .GT. SBIG) GO TO 13 ELSE IF ( NCD .EQ. 2 .AND. (LWK .LT. (ND+1)*NM1 .OR. . (PER .AND. LWK .LT. (ND+2)*NM1)) ) GO TO 12 SIG = 0.D0 ENDIF C C Compute the sequence of parameter values T. C SCURV = ND .EQ. 3 IF (.NOT. SCURV) THEN CALL ARCL2D (NN,X,Y, T,IERR) ELSE CALL ARCL3D (NN,X,Y,Z, T,IERR) ENDIF IF (IERR .GT. 0) GO TO 14 C C Initialize iteration count ITER, and store uniform C tension factors, or initialize SIGMA to zeros. C ITER = 0 DO 1 I = 1,NM1 SIGMA(I) = SIG 1 CONTINUE IF (NCD .EQ. 1) THEN C C NCD = 1. C IF (.NOT. PER) THEN CALL YPC1 (NN,T,X, XP,IERR) CALL YPC1 (NN,T,Y, YP,IERR) IF (SCURV) CALL YPC1 (NN,T,Z, ZP,IERR) ELSE CALL YPC1P (NN,T,X, XP,IERR) CALL YPC1P (NN,T,Y, YP,IERR) IF (SCURV) CALL YPC1P (NN,T,Z, ZP,IERR) ENDIF IF (IERR .NE. 0) GO TO 16 IF (.NOT. UNIFRM) THEN C C Call SIGS for UNIFRM = FALSE. C CALL SIGS (NN,T,X,XP,STOL, SIGMA, DSMAX,IERR) CALL SIGS (NN,T,Y,YP,STOL, SIGMA, DSMAX,IERR) IF (SCURV) CALL SIGS (NN,T,Z,ZP,STOL, SIGMA, . DSMAX,IERR) ENDIF GO TO 10 ENDIF C C NCD = 2. C IF (.NOT. PER) THEN C C Nonperiodic case: call YPC2 and test for IENDC invalid. C XP1 = XP(1) XPN = XP(NN) CALL YPC2 (NN,T,X,SIGMA,IENDC,IENDC,XP1,XPN, . WK, XP,IERR) YP1 = YP(1) YPN = YP(NN) CALL YPC2 (NN,T,Y,SIGMA,IENDC,IENDC,YP1,YPN, . WK, YP,IERR) IF (SCURV) THEN ZP1 = ZP(1) ZPN = ZP(NN) CALL YPC2 (NN,T,Z,SIGMA,IENDC,IENDC,ZP1,ZPN, . WK, ZP,IERR) ENDIF IF (IERR .EQ. 1) GO TO 11 IF (IERR .GT. 1) GO TO 16 ELSE C C Periodic fit: call YPC2P. C CALL YPC2P (NN,T,X,SIGMA,WK, XP,IERR) CALL YPC2P (NN,T,Y,SIGMA,WK, YP,IERR) IF (SCURV) CALL YPC2P (NN,T,Z,SIGMA,WK, ZP,IERR) IF (IERR .NE. 0) GO TO 16 ENDIF IF (UNIFRM) GO TO 10 C C Iterate on calls to SIGS and YPC2 (or YPC2P). The C first ND*(N-1) WK locations are used to store the C derivative estimates XP, YP, (and ZP) from the C previous iteration. IW1 is the first free WK location C following the stored derivatives. C C DYP is the maximum relative change in a component of XP, C YP, or ZP. C ICNT is the number of tension factors which were C increased by SIGS. C DSMAX is the maximum relative change in a component of C SIGMA. C IW1 = ND*NM1 + 1 DO 5 ITER = 1,MAXIT DYP = 0.D0 DO 2 I = 2,NM1 WK(I) = XP(I) WK(NM1+I) = YP(I) 2 CONTINUE IF (SCURV) THEN DO 3 I = 2,NM1 WK(N2M2+I) = ZP(I) 3 CONTINUE ENDIF CALL SIGS (NN,T,X,XP,STOL, SIGMA, DSMAX,ICNT) CALL SIGS (NN,T,Y,YP,STOL, SIGMA, DSMAX,ICNT) IF (SCURV) CALL SIGS (NN,T,Z,ZP,STOL, SIGMA, DSMAX, . ICNT) IF (.NOT. PER) THEN CALL YPC2 (NN,T,X,SIGMA,IENDC,IENDC,XP1,XPN, . WK(IW1), XP,IERR) CALL YPC2 (NN,T,Y,SIGMA,IENDC,IENDC,YP1,YPN, . WK(IW1), YP,IERR) IF (SCURV) CALL YPC2 (NN,T,Z,SIGMA,IENDC,IENDC, . ZP1,ZPN,WK(IW1), ZP,IERR) ELSE CALL YPC2P (NN,T,X,SIGMA,WK(IW1), XP,IERR) CALL YPC2P (NN,T,Y,SIGMA,WK(IW1), YP,IERR) IF (SCURV) CALL YPC2P (NN,T,Z,SIGMA,WK(IW1), ZP, . IERR) ENDIF EZ = 0.D0 DO 4 I = 2,NM1 EX = ABS(XP(I)-WK(I)) IF (WK(I) .NE. 0.D0) EX = EX/ABS(WK(I)) EY = ABS(YP(I)-WK(NM1+I)) IF (WK(NM1+I) .NE. 0.D0) EY = EY/ABS(WK(NM1+I)) IF (SCURV) THEN EZ = ABS(ZP(I)-WK(N2M2+I)) IF (WK(N2M2+I) .NE. 0.D0) . EZ = EZ/ABS(WK(N2M2+I)) ENDIF DYP = MAX(DYP,EX,EY,EZ) 4 CONTINUE IF (ICNT .EQ. 0 .OR. DYP .LE. DYPTOL) GO TO 6 5 CONTINUE ITER = MAXIT C C Store convergence parameters in WK. C 6 WK(1) = DYP WK(2) = DSMAX C C No error encountered. C 10 IER = ITER RETURN C C Invalid input parameter N, NCD, or IENDC. C 11 IER = -1 RETURN C C LWK too small. C 12 IER = -2 RETURN C C UNIFRM = TRUE and SIGMA(1) outside its valid range. C 13 IER = -3 RETURN C C Adjacent duplicate data points encountered. C 14 IER = -4 RETURN C C Error flag returned by YPC1, YPC1P, YPC2, or YPC2P: C T is not strictly increasing. C 16 IER = -6 RETURN END C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION TSINTL (A,B,N,X,Y,YP, . SIGMA, IER) INTEGER N, IER DOUBLE PRECISION A, B, X(N), Y(N), YP(N), SIGMA(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C This function computes the integral from A to B of a C Hermite interpolatory tension spline H. C C On input: C C A,B = Lower and upper limits of integration, re- C spectively. Note that -TSINTL(B,A,...) = C TSINTL(A,B,...). C C N = Number of data points. N .GE. 2. C C X = Array of length N containing the abscissae. C These must be in strictly increasing order: C X(I) < X(I+1) for I = 1,...,N-1. C C Y = Array of length N containing data values. C H(X(I)) = Y(I) for I = 1,...,N. C C YP = Array of length N containing first deriva- C tives. HP(X(I)) = YP(I) for I = 1,...,N, where C HP denotes the derivative of H. C C SIGMA = Array of length N-1 containing tension fac- C tors whose absolute values determine the C balance between cubic and linear in each C interval. SIGMA(I) is associated with int- C erval (I,I+1) for I = 1,...,N-1. C C Input parameters are not altered by this function. C C On output: C C IER = Error indicator: C IER = 0 if no errors were encountered and C X(1) .LE. T .LE. X(N) for T = A and C T = B, or A = B. C IER = 1 if no errors were encountered but C extrapolation was necessary: A or B C not in the interval (X(1),X(N)). C IER = -1 IF N < 2. C IER = -2 if the abscissae are not in strictly C increasing order. Only those in or C adjacent to the interval of integra- C tion are tested. C C TSINTL = Integral of H from A to B, or zero if C IER < 0. C C Modules required by TSINTL: INTRVL, SNHCSH C C Intrinsic functions called by TSINTL: ABS, EXP, MAX, MIN C C*********************************************************** C INTEGER I, IL, ILP1, IMAX, IMIN, IP1, IU, IUP1 DOUBLE PRECISION B1, B2, CM, CM1, CM2, CMM, CMM1, . CMM2, D1, D2, DX, E, E1, E2, EMS, S, . S1, S2, SB1, SB2, SBIG, SIG, SM, SM1, . SM2, SUM, T, TM, TP, U, XL, XU, Y1, . Y2 INTEGER INTRVL C DATA SBIG/85.D0/ IF (N .LT. 2) GO TO 7 C C Accumulate the integral from XL to XU in SUM. C XL = MIN(A,B) XU = MAX(A,B) SUM = 0.D0 IER = 0 IF (XL .EQ. XU) GO TO 6 C C Find left-end indexes of intervals containing XL and XU. C If XL < X(1) or XU > X(N), extrapolation is performed C using the leftmost or rightmost interval. C IL = INTRVL (XL,N,X) IU = INTRVL (XU,N,X) IF (XL .LT. X(1) .OR. XU .GT. X(N)) IER = 1 ILP1 = IL + 1 IMIN = IL IF (XL .EQ. X(IL)) GO TO 2 C C Compute the integral from XL to X(IL+1). C DX = X(ILP1) - X(IL) IF (DX .LE. 0.D0) GO TO 8 U = X(ILP1) - XL IF (U .EQ. 0.D0) GO TO 1 B1 = U/DX Y2 = Y(ILP1) S = (Y2-Y(IL))/DX S2 = YP(ILP1) D1 = S - YP(IL) D2 = S2 - S SIG = ABS(SIGMA(IL)) IF (SIG .LT. 1.D-9) THEN C C SIG = 0. C SUM = SUM + U*(Y2 - U*(6.D0*S2 - B1*(4.D0*D2 + . (3.D0*B1-4.D0)*(D1-D2)))/12.D0) ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5. C SB1 = SIG*B1 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB1, SM1,CM1,CMM1) E = SIG*SM - CMM - CMM SUM = SUM + U*(Y2 - S2*U/2.D0) + ((CM*CMM1-SM*SM1)* . (D1+D2) + SIG*(CM*SM1-(SM+SIG)*CMM1)*D2)/ . ((SIG/DX)**2*E) ELSE C C SIG > .5. C SB1 = SIG*B1 SB2 = SIG - SB1 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN SUM = SUM + U*(Y2 - S*U/2.D0) ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1.D0 - EMS TP = 1.D0 + EMS T = SB1*SB1/2.D0 + 1.D0 E = TM*(SIG*TP - TM - TM) SUM = SUM +U*(Y2 - S2*U/2.D0)+(SIG*TM*(TP*T-E1-E2- . TM*SB1)*D2 - (TM*(TM*T-E1+E2-TP*SB1) + . SIG*(E1*EMS-E2+2.D0*SB1*EMS))*(D1+D2))/ . ((SIG/DX)**2*E) ENDIF ENDIF C C Test for XL and XU in the same interval. C 1 IMIN = ILP1 IF (IL .EQ. IU) GO TO 5 C C Add in the integral from X(IMIN) to X(J) for J = C Max(IL+1,IU). C 2 IMAX = MAX(IL,IU-1) DO 3 I = IMIN,IMAX IP1 = I + 1 DX = X(IP1) - X(I) IF (DX .LE. 0.D0) GO TO 8 SIG = ABS(SIGMA(I)) IF (SIG .LT. 1.D-9) THEN C C SIG = 0. C SUM = SUM + DX*((Y(I)+Y(IP1))/2.D0 - . DX*(YP(IP1)-YP(I))/12.D0) ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5. C CALL SNHCSH (SIG, SM,CM,CMM) E = SIG*SM - CMM - CMM SUM = SUM + DX*(Y(I)+Y(IP1) - DX*E*(YP(IP1)-YP(I)) . /(SIG*SIG*CM))/2.D0 ELSE C C SIG > .5. C EMS = EXP(-SIG) SUM = SUM + DX*(Y(I)+Y(IP1) - DX*(SIG*(1.D0+EMS)/ . (1.D0-EMS)-2.D0)*(YP(IP1)-YP(I))/ . (SIG*SIG))/2.D0 ENDIF 3 CONTINUE C C Add in the integral from X(IU) to XU if IU > IL. C IF (IL .EQ. IU) GO TO 4 IUP1 = IU + 1 DX = X(IUP1) - X(IU) IF (DX .LE. 0.D0) GO TO 8 U = XU - X(IU) IF (U .EQ. 0.D0) GO TO 6 B2 = U/DX Y1 = Y(IU) S = (Y(IUP1)-Y1)/DX S1 = YP(IU) D1 = S - S1 D2 = YP(IUP1) - S SIG = ABS(SIGMA(IU)) IF (SIG .LT. 1.D-9) THEN C C SIG = 0. C SUM = SUM + U*(Y1 + U*(6.D0*S1 + B2*(4.D0*D1 + . (4.D0-3.D0*B2)*(D1-D2)))/12.D0) ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5. C SB2 = SIG*B2 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB2, SM2,CM2,CMM2) E = SIG*SM - CMM - CMM SUM = SUM + U*(Y1 + S1*U/2.D0) + ((CM*CMM2-SM*SM2)* . (D1+D2) + SIG*(CM*SM2-(SM+SIG)*CMM2)*D1) . /((SIG/DX)**2*E) ELSE C C SIG > .5. C SB2 = SIG*B2 SB1 = SIG - SB2 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN SUM = SUM + U*(Y1 + S*U/2.D0) ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1.D0 - EMS TP = 1.D0 + EMS T = SB2*SB2/2.D0 + 1.D0 E = TM*(SIG*TP - TM - TM) SUM = SUM +U*(Y1 + S1*U/2.D0)+(SIG*TM*(TP*T-E1-E2- . TM*SB2)*D1 - (TM*(TM*T-E2+E1-TP*SB2) + . SIG*(E2*EMS-E1+2.D0*SB2*EMS))*(D1+D2))/ . ((SIG/DX)**2*E) ENDIF ENDIF GO TO 6 C C IL = IU and SUM contains the integral from XL to X(IL+1). C Subtract off the integral from XU to X(IL+1). DX and C SIG were computed above. C 4 Y2 = Y(ILP1) S = (Y2-Y(IL))/DX S2 = YP(ILP1) D1 = S - YP(IL) D2 = S2 - S C 5 U = X(ILP1) - XU IF (U .EQ. 0.D0) GO TO 6 B1 = U/DX IF (SIG .LT. 1.D-9) THEN C C SIG = 0. C SUM = SUM - U*(Y2 - U*(6.D0*S2 - B1*(4.D0*D2 + . (3.D0*B1-4.D0)*(D1-D2)))/12.D0) ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5. C SB1 = SIG*B1 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB1, SM1,CM1,CMM1) E = SIG*SM - CMM - CMM SUM = SUM - U*(Y2 - S2*U/2.D0) - ((CM*CMM1-SM*SM1)* . (D1+D2) + SIG*(CM*SM1-(SM+SIG)*CMM1)*D2) . /((SIG/DX)**2*E) ELSE C C SIG > .5. C SB1 = SIG*B1 SB2 = SIG - SB1 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN SUM = SUM - U*(Y2 - S*U/2.D0) ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1.D0 - EMS TP = 1.D0 + EMS T = SB1*SB1/2.D0 + 1.D0 E = TM*(SIG*TP - TM - TM) SUM = SUM -U*(Y2 - S2*U/2.D0)-(SIG*TM*(TP*T-E1-E2- . TM*SB1)*D2 - (TM*(TM*T-E1+E2-TP*SB1) + . SIG*(E1*EMS-E2+2.D0*SB1*EMS))*(D1+D2))/ . ((SIG/DX)**2*E) ENDIF ENDIF C C No errors were encountered. Adjust the sign of SUM. C 6 IF (XL .EQ. B) SUM = -SUM TSINTL = SUM RETURN C C N < 2. C 7 IER = -1 TSINTL = 0.D0 RETURN C C Abscissae not strictly increasing. C 8 IER = -2 TSINTL = 0.D0 RETURN END C ---------------------------------------------------------- SUBROUTINE ARCL2D (N,X,Y, T,IER) INTEGER N, IER DOUBLE PRECISION X(N), Y(N), T(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C Given an ordered sequence of points (X,Y) defining a C polygonal curve in the plane, this subroutine computes the C sequence T of cumulative arc lengths along the curve: C T(1) = 0 and, for 2 .LE. K .LE. N, T(K) is the sum of C Euclidean distances between (X(I-1),Y(I-1)) and (X(I),Y(I)) C for I = 2,...,K. A closed curve corresponds to X(1) = C X(N) and Y(1) = Y(N), and more generally, duplicate points C are permitted but must not be adjacent. Thus, T contains C a strictly increasing sequence of values which may be used C as parameters for fitting a smooth curve to the sequence C of points. C C On input: C C N = Number of points defining the curve. N .GE. 2. C C X,Y = Arrays of length N containing the coordinates C of the points. C C The above parameters are not altered by this routine. C C T = Array of length at least N. C C On output: C C T = Array containing cumulative arc lengths defined C above unless IER > 0. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N < 2. C IER = I if X(I-1) = X(I) and Y(I-1) = Y(I) for C some I in the range 2,...,N. C C Modules required by ARCL2D: None C C Intrinsic function called by ARCL2D: SQRT C C*********************************************************** C INTEGER I, NN DOUBLE PRECISION DS C NN = N IF (NN .LT. 2) GO TO 2 T(1) = 0.D0 DO 1 I = 2,NN DS = (X(I)-X(I-1))**2 + (Y(I)-Y(I-1))**2 IF (DS .EQ. 0.D0) GO TO 3 T(I) = T(I-1) + SQRT(DS) 1 CONTINUE IER = 0 RETURN C C N is outside its valid range. C 2 IER = 1 RETURN C C Points I-1 and I coincide. C 3 IER = I RETURN END C ---------------------------------------------------------- SUBROUTINE ARCL3D (N,X,Y,Z, T,IER) INTEGER N, IER DOUBLE PRECISION X(N), Y(N), Z(N), T(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C Given an ordered sequence of points (X,Y,Z) defining a C polygonal curve in 3-space, this subroutine computes the C sequence T of cumulative arc lengths along the curve: C T(1) = 0 and, for 2 .LE. K .LE. N, T(K) is the sum of C Euclidean distances between (X(I-1),Y(I-1),Z(I-1)) and C (X(I),Y(I),Z(I)) for I = 2,...,K. A closed curve corre- C sponds to X(1) = X(N), Y(1) = Y(N), and Z(1) = Z(N). More C generally, duplicate points are permitted but must not be C adjacent. Thus, T contains a strictly increasing sequence C of values which may be used as parameters for fitting a C smooth curve to the sequence of points. C C On input: C C N = Number of points defining the curve. N .GE. 2. C C X,Y,Z = Arrays of length N containing the coordi- C nates of the points. C C The above parameters are not altered by this routine. C C T = Array of length at least N. C C On output: C C T = Array containing cumulative arc lengths defined C above unless IER > 0. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N < 2. C IER = I if X(I-1) = X(I), Y(I-1) = Y(I), and C Z(I-1) = Z(I) for some I in the range C 2,...,N. C C Modules required by ARCL3D: None C C Intrinsic function called by ARCL3D: SQRT C C*********************************************************** C INTEGER I, NN DOUBLE PRECISION DS C NN = N IF (NN .LT. 2) GO TO 2 T(1) = 0.D0 DO 1 I = 2,NN DS = (X(I)-X(I-1))**2 + (Y(I)-Y(I-1))**2 + . (Z(I)-Z(I-1))**2 IF (DS .EQ. 0.D0) GO TO 3 T(I) = T(I-1) + SQRT(DS) 1 CONTINUE IER = 0 RETURN C C N is outside its valid range. C 2 IER = 1 RETURN C C Points I-1 and I coincide. C 3 IER = I RETURN END C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION ENDSLP (X1,X2,X3,Y1,Y2,Y3, . SIGMA) DOUBLE PRECISION X1, X2, X3, Y1, Y2, Y3, SIGMA C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C Given data values associated with a strictly increasing C or decreasing sequence of three abscissae X1, X2, and X3, C this function returns a derivative estimate at X1 based on C the tension spline H(x) which interpolates the data points C and has third derivative equal to zero at X1. Letting S1 C denote the slope defined by the first two points, the est- C mate is obtained by constraining the derivative of H at X1 C so that it has the same sign as S1 and its magnitude is C at most 3*abs(S1). If SIGMA = 0, H(x) is quadratic and C the derivative estimate is identical to the value computed C by Subroutine YPC1 at the first point (or the last point C if the abscissae are decreasing). C C On input: C C X1,X2,X3 = Abscissae satisfying either X1 < X2 < X3 C or X1 > X2 > X3. C C Y1,Y2,Y3 = Data values associated with the abscis- C sae. H(X1) = Y1, H(X2) = Y2, and H(X3) C = Y3. C C SIGMA = Tension factor associated with H in inter- C val (X1,X2) or (X2,X1). C C Input parameters are not altered by this function. C C On output: C C ENDSLP = (Constrained) derivative of H at X1, or C zero if the abscissae are not strictly C monotonic. C C Module required by ENDSLP: SNHCSH C C Intrinsic functions called by ENDSLP: ABS, EXP, MAX, MIN C C*********************************************************** C DOUBLE PRECISION COSHM1, COSHMS, DUMMY, DX1, DXS, S1, . SIG1, SIGS, T C DX1 = X2 - X1 DXS = X3 - X1 IF (DX1*(DXS-DX1) .LE. 0.D0) GO TO 2 SIG1 = ABS(SIGMA) IF (SIG1 .LT. 1.D-9) THEN C C SIGMA = 0: H is the quadratic interpolant. C T = (DX1/DXS)**2 GO TO 1 ENDIF SIGS = SIG1*DXS/DX1 IF (SIGS .LE. .5D0) THEN C C 0 < SIG1 < SIGS .LE. .5: compute approximations to C COSHM1 = COSH(SIG1)-1 and COSHMS = COSH(SIGS)-1. C CALL SNHCSH (SIG1, DUMMY,COSHM1,DUMMY) CALL SNHCSH (SIGS, DUMMY,COSHMS,DUMMY) T = COSHM1/COSHMS ELSE C C SIGS > .5: compute T = COSHM1/COSHMS. C T = EXP(SIG1-SIGS)*((1.D0-EXP(-SIG1))/ . (1.D0-EXP(-SIGS)))**2 ENDIF C C The derivative of H at X1 is C T = ((Y3-Y1)*COSHM1-(Y2-Y1)*COSHMS)/ C (DXS*COSHM1-DX1*COSHMS). C C ENDSLP = T unless T*S1 < 0 or abs(T) > 3*abs(S1). C 1 T = ((Y3-Y1)*T-Y2+Y1)/(DXS*T-DX1) S1 = (Y2-Y1)/DX1 IF (S1 .GE. 0.D0) THEN ENDSLP = MIN(MAX(0.D0,T), 3.D0*S1) ELSE ENDSLP = MAX(MIN(0.D0,T), 3.D0*S1) ENDIF RETURN C C Error in the abscissae. C 2 ENDSLP = 0.D0 RETURN END C ---------------------------------------------------------- SUBROUTINE SIGS (N,X,Y,YP,TOL, SIGMA, DSMAX,IER) INTEGER N, IER DOUBLE PRECISION X(N), Y(N), YP(N), TOL, SIGMA(N), . DSMAX C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C Given a set of abscissae X with associated data values Y C and derivatives YP, this subroutine determines the small- C est (nonnegative) tension factors SIGMA such that the Her- C mite interpolatory tension spline H(x) preserves local C shape properties of the data. In an interval (X1,X2) with C data values Y1,Y2 and derivatives YP1,YP2, the properties C of the data are C C Monotonicity: S, YP1, and YP2 are nonnegative or C nonpositive, C and C Convexity: YP1 .LE. S .LE. YP2 or YP1 .GE. S C .GE. YP2, C C where S = (Y2-Y1)/(X2-X1). The corresponding properties C of H are constant sign of the first and second deriva- C tives, respectively. Note that, unless YP1 = S = YP2, in- C finite tension is required (and H is linear on the inter- C val) if S = 0 in the case of monotonicity, or if YP1 = S C or YP2 = S in the case of convexity. C C SIGS may be used in conjunction with Subroutine YPC2 C (or YPC2P) in order to produce a C-2 interpolant which C preserves the shape properties of the data. This is C achieved by calling YPC2 with SIGMA initialized to the C zero vector, and then alternating calls to SIGS with C calls to YPC2 until the change in SIGMA is small (refer to C the parameter descriptions for SIGMA, DSMAX and IER), or C the maximum relative change in YP is bounded by a toler- C ance (a reasonable value is .01). A similar procedure may C be used to produce a C-2 shape-preserving smoothing curve C (Subroutine SMCRV). C C Refer to Subroutine SIGBI for a means of selecting mini- C mum tension factors to satisfy more general constraints. C C On input: C C N = Number of data points. N .GE. 2. C C X = Array of length N containing a strictly in- C creasing sequence of abscissae: X(I) < X(I+1) C for I = 1,...,N-1. C C Y = Array of length N containing data values (or C function values computed by SMCRV) associated C with the abscissae. H(X(I)) = Y(I) for I = C 1,...,N. C C YP = Array of length N containing first derivatives C of H at the abscissae. Refer to Subroutines C YPC1, YPC1P, YPC2, YPC2P, and SMCRV. C C TOL = Tolerance whose magnitude determines how close C each tension factor is to its optimal value C when nonzero finite tension is necessary and C sufficient to satisfy the constraint: C abs(TOL) is an upper bound on the magnitude C of the smallest (nonnegative) or largest (non- C positive) value of the first or second deriva- C tive of H in the interval. Thus, the con- C straint is satisfied, but possibly with more C tension than necessary. TOL should be set to C 0 for optimal tension. C C The above parameters are not altered by this routine. C C SIGMA = Array of length N-1 containing minimum val- C ues of the tension factors. SIGMA(I) is as- C sociated with interval (I,I+1) and SIGMA(I) C .GE. 0 for I = 1,...,N-1. SIGMA should be C set to the zero vector if minimal tension C is desired, and should be unchanged from a C previous call in order to ensure convergence C of the C-2 iterative procedure. C C On output: C C SIGMA = Array containing tension factors for which C H(x) preserves the properties of the data, C with the restriction that SIGMA(I) .LE. 85 C for all I (unless the input value is larger). C The factors are as small as possible (within C the tolerance), but not less than their C input values. If infinite tension is re- C quired in interval (X(I),X(I+1)), then C SIGMA(I) = 85 (and H is an approximation to C the linear interpolant on the interval), C and if neither property is satisfied by the C data, then SIGMA(I) = 0 (unless the input C value is positive), and thus H is cubic in C the interval. C C DSMAX = Maximum increase in a component of SIGMA C from its input value. The increase is a C relative change if the input value is C nonzero, and an absolute change otherwise. C C IER = Error indicator and information flag: C IER = I if no errors were encountered and I C components of SIGMA were altered from C their input values for 0 .LE. I .LE. C N-1. C IER = -1 if N < 2. SIGMA is not altered in C this case. C IER = -I if X(I) .LE. X(I-1) for some I in the C range 2,...,N. SIGMA(J-1) is unal- C tered for J = I,...,N in this case. C C Modules required by SIGS: SNHCSH, STORE C C Intrinsic functions called by SIGS: ABS, EXP, MAX, MIN, C SIGN, SQRT C C*********************************************************** C INTEGER I, ICNT, IP1, LUN, NIT, NM1, NITSTEP, NITMAX DOUBLE PRECISION A, C1, C2, COSHM, COSHMM, D0, D1, . D1D2, D1PD2, D2, DMAX, DSIG, DSM, DX, . E, EMS, EMS2, F, F0, FMAX, FNEG, FP, . FTOL, RTOL, S, S1, S2, SBIG, SCM, . SGN, SIG, SIGIN, SINHM, SSINH, SSM, . STOL, T, T0, T1, T2, TM, TP1 DOUBLE PRECISION STORE C C The following line altered by Pat Scott Jan 31 2008 DATA SBIG/85.D0/, LUN/-1/, NITSTEP/200/, NITMAX/500/ NM1 = N - 1 IF (NM1 .LT. 1) GO TO 9 C C Compute an absolute tolerance FTOL = abs(TOL) and a C relative tolerance RTOL = 100*MACHEPS (now 200*MACHEPS). C FTOL = ABS(TOL) RTOL = 1.D0 1 RTOL = RTOL/2.D0 IF (STORE(RTOL+1.D0) .GT. 1.D0) GO TO 1 C The following line altered by Pat Scott Jan 31 2008 RTOL = RTOL*400.D0 C C Initialize change counter ICNT and maximum change DSM for C loop on intervals. C ICNT = 0 DSM = 0.D0 DO 8 I = 1,NM1 IF (LUN .GE. 0) WRITE (LUN,100) I 100 FORMAT (//1X,'SIGS -- INTERVAL',I4) IP1 = I + 1 DX = X(IP1) - X(I) IF (DX .LE. 0.D0) GO TO 10 SIGIN = SIGMA(I) IF (SIGIN .GE. SBIG) GO TO 8 C C Compute first and second differences. C S1 = YP(I) S2 = YP(IP1) S = (Y(IP1)-Y(I))/DX D1 = S - S1 D2 = S2 - S D1D2 = D1*D2 C C Test for infinite tension required to satisfy either C property. C SIG = SBIG IF ((D1D2 .EQ. 0.D0 .AND. S1 .NE. S2) .OR. . (S .EQ. 0.D0 .AND. S1*S2 .GT. 0.D0)) GO TO 7 C C Test for SIGMA = 0 sufficient. The data satisfies convex- C ity iff D1D2 .GE. 0, and D1D2 = 0 implies S1 = S = S2. C SIG = 0.D0 IF (D1D2 .LT. 0.D0) GO TO 3 IF (D1D2 .EQ. 0.D0) GO TO 7 T = MAX(D1/D2,D2/D1) IF (T .LE. 2.D0) GO TO 7 TP1 = T + 1.D0 C C Convexity: Find a zero of F(SIG) = SIG*COSHM(SIG)/ C SINHM(SIG) - TP1. C C F(0) = 2-T < 0, F(TP1) .GE. 0, the derivative of F C vanishes at SIG = 0, and the second derivative of F is C .2 at SIG = 0. A quadratic approximation is used to C obtain a starting point for the Newton method. C SIG = SQRT(10.D0*T-20.D0) NIT = 0 C C Top of loop: C 2 IF (SIG .LE. .5D0) THEN CALL SNHCSH (SIG, SINHM,COSHM,COSHMM) T1 = COSHM/SINHM FP = T1 + SIG*(SIG/SINHM - T1*T1 + 1.D0) ELSE C C Scale SINHM and COSHM by 2*EXP(-SIG) in order to avoid C overflow with large SIG. C EMS = EXP(-SIG) SSM = 1.D0 - EMS*(EMS+SIG+SIG) T1 = (1.D0-EMS)*(1.D0-EMS)/SSM FP = T1 + SIG*(2.D0*SIG*EMS/SSM - T1*T1 + 1.D0) ENDIF C F = SIG*T1 - TP1 IF (LUN .GE. 0) WRITE (LUN,110) SIG, F, FP 110 FORMAT (5X,'CONVEXITY -- SIG = ',D15.8, . ', F(SIG) = ',D15.8/1X,35X,'FP(SIG) = ', . D15.8) NIT = NIT + 1 C C Test for convergence. C IF (FP .LE. 0.D0) GO TO 7 DSIG = -F/FP IF (ABS(DSIG) .LE. RTOL*SIG .OR. (F .GE. 0.D0 . .AND. F .LE. FTOL) .OR. ABS(F) .LE. RTOL) . GO TO 7 C C Update SIG. C SIG = SIG + DSIG C The following lines added by Pat Scott Jan 31 2008 IF (NIT .EQ. NITSTEP) THEN RTOL = RTOL * 2.D0 WRITE(*,*) 'WARNING: RTOL in SIGS doubled to aid convergence' ENDIF IF (NIT .EQ. NITMAX) THEN WRITE(*,*) 'Error: SIGS did not converge' GO TO 9 ENDIF C end addition on Jan 31 2008 GO TO 2 C C Convexity cannot be satisfied. Monotonicity can be satis- C fied iff S1*S .GE. 0 and S2*S .GE. 0 since S .NE. 0. C 3 IF (S1*S .LT. 0.D0 .OR. S2*S .LT. 0.D0) GO TO 7 T0 = 3.D0*S - S1 - S2 D0 = T0*T0 - S1*S2 C C SIGMA = 0 is sufficient for monotonicity iff S*T0 .GE. 0 C or D0 .LE. 0. C IF (D0 .LE. 0.D0 .OR. S*T0 .GE. 0.D0) GO TO 7 C C Monotonicity: find a zero of F(SIG) = SIGN(S)*HP(R), C where HPP(R) = 0 and HP, HPP denote derivatives of H. C F has a unique zero, F(0) < 0, and F approaches abs(S) C as SIG increases. C C Initialize parameters for the secant method. The method C uses three points: (SG0,F0), (SIG,F), and C (SNEG,FNEG), where SG0 and SNEG are defined implicitly C by DSIG = SIG - SG0 and DMAX = SIG - SNEG. C SGN = SIGN(1.D0,S) SIG = SBIG FMAX = SGN*(SIG*S-S1-S2)/(SIG-2.D0) IF (FMAX .LE. 0.D0) GO TO 7 STOL = RTOL*SIG F = FMAX F0 = SGN*D0/(3.D0*(D1-D2)) FNEG = F0 DSIG = SIG DMAX = SIG D1PD2 = D1 + D2 NIT = 0 C C Top of loop: compute the change in SIG by linear C interpolation. C 4 DSIG = -F*DSIG/(F-F0) IF (LUN .GE. 0) WRITE (LUN,120) DSIG 120 FORMAT (5X,'MONOTONICITY -- DSIG = ',D15.8) IF ( ABS(DSIG) .GT. ABS(DMAX) .OR. . DSIG*DMAX .GT. 0. ) GO TO 6 C C Restrict the step-size such that abs(DSIG) .GE. STOL/2. C Note that DSIG and DMAX have opposite signs. C IF (ABS(DSIG) .LT. STOL/2.D0) . DSIG = -SIGN(STOL/2.D0,DMAX) C C Update SIG, F0, and F. C SIG = SIG + DSIG F0 = F IF (SIG .LE. .5D0) THEN C C Use approximations to the hyperbolic functions designed C to avoid cancellation error with small SIG. C CALL SNHCSH (SIG, SINHM,COSHM,COSHMM) C1 = SIG*COSHM*D2 - SINHM*D1PD2 C2 = SIG*(SINHM+SIG)*D2 - COSHM*D1PD2 A = C2 - C1 E = SIG*SINHM - COSHMM - COSHMM ELSE C C Scale SINHM and COSHM by 2*EXP(-SIG) in order to avoid C overflow with large SIG. C EMS = EXP(-SIG) EMS2 = EMS + EMS TM = 1.D0 - EMS SSINH = TM*(1.D0+EMS) SSM = SSINH - SIG*EMS2 SCM = TM*TM C1 = SIG*SCM*D2 - SSM*D1PD2 C2 = SIG*SSINH*D2 - SCM*D1PD2 C C R is in (0,1) and well-defined iff HPP(X1)*HPP(X2) < 0. C F = FMAX IF (C1*(SIG*SCM*D1 - SSM*D1PD2) .GE. 0.D0) GO TO 5 A = EMS2*(SIG*TM*D2 + (TM-SIG)*D1PD2) IF (A*(C2+C1) .LT. 0.D0) GO TO 5 E = SIG*SSINH - SCM - SCM ENDIF C F = (SGN*(E*S2-C2) + SQRT(A*(C2+C1)))/E C C Update number of iterations NIT. C 5 NIT = NIT + 1 IF (LUN .GE. 0) WRITE (LUN,130) NIT, SIG, F 130 FORMAT (1X,10X,I2,' -- SIG = ',D15.8,', F = ', . D15.8) C C Test for convergence. C STOL = RTOL*SIG IF ( ABS(DMAX) .LE. STOL .OR. (F .GE. 0.D0 .AND. . F .LE. FTOL) .OR. ABS(F) .LE. RTOL ) GO TO 7 DMAX = DMAX + DSIG IF ( F0*F .GT. 0.D0 .AND. ABS(F) .GE. ABS(F0) ) . GO TO 6 IF (F0*F .LE. 0.D0) THEN C C F and F0 have opposite signs. Update (SNEG,FNEG) to C (SG0,F0) so that F and FNEG always have opposite C signs. If SIG is closer to SNEG than SG0 and abs(F) < C abs(FNEG), then swap (SNEG,FNEG) with (SG0,F0). C T1 = DMAX T2 = FNEG DMAX = DSIG FNEG = F0 IF ( ABS(DSIG) .GT. ABS(T1) .AND. . ABS(F) .LT. ABS(T2) ) THEN C DSIG = T1 F0 = T2 ENDIF ENDIF GO TO 4 C C Bottom of loop: F0*F > 0 and the new estimate would C be outside of the bracketing interval of length C abs(DMAX). Reset (SG0,F0) to (SNEG,FNEG). C 6 DSIG = DMAX F0 = FNEG GO TO 4 C C Update SIGMA(I), ICNT, and DSM if necessary. C 7 SIG = MIN(SIG,SBIG) IF (SIG .GT. SIGIN) THEN SIGMA(I) = SIG ICNT = ICNT + 1 DSIG = SIG-SIGIN IF (SIGIN .GT. 0.D0) DSIG = DSIG/SIGIN DSM = MAX(DSM,DSIG) ENDIF 8 CONTINUE C C No errors encountered. C DSMAX = DSM IER = ICNT RETURN C C N < 2. C 9 DSMAX = 0.D0 IER = -1 RETURN C C X(I+1) .LE. X(I). C 10 DSMAX = DSM IER = -IP1 RETURN END C ---------------------------------------------------------- SUBROUTINE SNHCSH (X, SINHM,COSHM,COSHMM) DOUBLE PRECISION X, SINHM, COSHM, COSHMM C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/20/96 C C This subroutine computes approximations to the modified C hyperbolic functions defined below with relative error C bounded by 3.4E-20 for a floating point number system with C sufficient precision. C C Note that the 21-digit constants in the data statements C below may not be acceptable to all compilers. C C On input: C C X = Point at which the functions are to be C evaluated. C C X is not altered by this routine. C C On output: C C SINHM = sinh(X) - X. C C COSHM = cosh(X) - 1. C C COSHMM = cosh(X) - 1 - X*X/2. C C Modules required by SNHCSH: None C C Intrinsic functions called by SNHCSH: ABS, EXP C C*********************************************************** C DOUBLE PRECISION AX, EXPX, F, P, P1, P2, P3, P4, Q, . Q1, Q2, Q3, Q4, XC, XS, XSD2, XSD4 C DATA P1/-3.51754964808151394800D5/, . P2/-1.15614435765005216044D4/, . P3/-1.63725857525983828727D2/, . P4/-7.89474443963537015605D-1/ DATA Q1/-2.11052978884890840399D6/, . Q2/3.61578279834431989373D4/, . Q3/-2.77711081420602794433D2/, . Q4/1.D0/ AX = ABS(X) XS = AX*AX IF (AX .LE. .5D0) THEN C C Approximations for small X: C XC = X*XS P = ((P4*XS+P3)*XS+P2)*XS+P1 Q = ((Q4*XS+Q3)*XS+Q2)*XS+Q1 SINHM = XC*(P/Q) XSD4 = .25D0*XS XSD2 = XSD4 + XSD4 P = ((P4*XSD4+P3)*XSD4+P2)*XSD4+P1 Q = ((Q4*XSD4+Q3)*XSD4+Q2)*XSD4+Q1 F = XSD4*(P/Q) COSHMM = XSD2*F*(F+2.D0) COSHM = COSHMM + XSD2 ELSE C C Approximations for large X: C EXPX = EXP(AX) SINHM = -(((1.D0/EXPX+AX)+AX)-EXPX)/2.D0 IF (X .LT. 0.D0) SINHM = -SINHM COSHM = ((1.D0/EXPX-2.D0)+EXPX)/2.D0 COSHMM = COSHM - XS/2.D0 ENDIF RETURN END C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION STORE (X) DOUBLE PRECISION X C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This function forces its argument X to be stored in a C memory location, thus providing a means of determining C floating point number characteristics (such as the machine C precision) when it is necessary to avoid computation in C high precision registers. C C On input: C C X = Value to be stored. C C X is not altered by this function. C C On output: C C STORE = Value of X after it has been stored and C possibly truncated or rounded to the single C precision word length. C C Modules required by STORE: None C C*********************************************************** C DOUBLE PRECISION Y COMMON/STCOM/Y Y = X STORE = Y RETURN END C ---------------------------------------------------------- SUBROUTINE YPCOEF (SIGMA,DX, D,SD) DOUBLE PRECISION SIGMA, DX, D, SD C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C This subroutine computes the coefficients of the deriva- C tives in the symmetric diagonally dominant tridiagonal C system associated with the C-2 derivative estimation pro- C cedure for a Hermite interpolatory tension spline. C C On input: C C SIGMA = Nonnegative tension factor associated with C an interval. C C DX = Positive interval width. C C Input parameters are not altered by this routine. C C On output: C C D = Component of the diagonal term associated with C the interval. D = SIG*(SIG*COSHM(SIG) - C SINHM(SIG))/(DX*E), where SIG = SIGMA and E = C SIG*SINH(SIG) - 2*COSHM(SIG). C C SD = Subdiagonal (superdiagonal) term. SD = SIG* C SINHM(SIG)/E. C C Module required by YPCOEF: SNHCSH C C Intrinsic function called by YPCOEF: EXP C C*********************************************************** C DOUBLE PRECISION COSHM, COSHMM, E, EMS, SCM, SIG, . SINHM, SSINH, SSM C SIG = SIGMA IF (SIG .LT. 1.D-9) THEN C C SIG = 0: cubic interpolant. C D = 4.D0/DX SD = 2.D0/DX ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5: use approximations designed to avoid C cancellation error in the hyperbolic C functions when SIGMA is small. C CALL SNHCSH (SIG, SINHM,COSHM,COSHMM) E = (SIG*SINHM - COSHMM - COSHMM)*DX D = SIG*(SIG*COSHM-SINHM)/E SD = SIG*SINHM/E ELSE C C SIG > .5: scale SINHM and COSHM by 2*EXP(-SIG) in order C to avoid overflow when SIGMA is large. C EMS = EXP(-SIG) SSINH = 1.D0 - EMS*EMS SSM = SSINH - 2.D0*SIG*EMS SCM = (1.D0-EMS)*(1.D0-EMS) E = (SIG*SSINH - SCM - SCM)*DX D = SIG*(SIG*SCM-SSM)/E SD = SIG*SSM/E ENDIF RETURN END C ---------------------------------------------------------- SUBROUTINE YPC1 (N,X,Y, YP,IER) INTEGER N, IER DOUBLE PRECISION X(N), Y(N), YP(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This subroutine employs a three-point quadratic interpo- C lation method to compute local derivative estimates YP C associated with a set of data points. The interpolation C formula is the monotonicity-constrained parabolic method C described in the reference cited below. A Hermite int- C erpolant of the data values and derivative estimates pre- C serves monotonicity of the data. Linear interpolation is C used if N = 2. The method is invariant under a linear C scaling of the coordinates but is not additive. C C On input: C C N = Number of data points. N .GE. 2. C C X = Array of length N containing a strictly in- C creasing sequence of abscissae: X(I) < X(I+1) C for I = 1,...,N-1. C C Y = Array of length N containing data values asso- C ciated with the abscissae. C C Input parameters are not altered by this routine. C C On output: C C YP = Array of length N containing estimated deriv- C atives at the abscissae unless IER .NE. 0. C YP is not altered if IER = 1, and is only par- C tially defined if IER > 1. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N < 2. C IER = I if X(I) .LE. X(I-1) for some I in the C range 2,...,N. C C Reference: J. M. Hyman, "Accurate Monotonicity-preserving C Cubic Interpolation", LA-8796-MS, Los C Alamos National Lab, Feb. 1982. C C Modules required by YPC1: None C C Intrinsic functions called by YPC1: ABS, MAX, MIN, SIGN C C*********************************************************** C INTEGER I, NM1 DOUBLE PRECISION ASI, ASIM1, DX2, DXI, DXIM1, S2, SGN, . SI, SIM1, T C NM1 = N - 1 IF (NM1 .LT. 1) GO TO 2 I = 1 DXI = X(2) - X(1) IF (DXI .LE. 0.D0) GO TO 3 SI = (Y(2)-Y(1))/DXI IF (NM1 .EQ. 1) THEN C C Use linear interpolation for N = 2. C YP(1) = SI YP(2) = SI IER = 0 RETURN ENDIF C C N .GE. 3. YP(1) = S1 + DX1*(S1-S2)/(DX1+DX2) unless this C results in YP(1)*S1 .LE. 0 or abs(YP(1)) > 3*abs(S1). C I = 2 DX2 = X(3) - X(2) IF (DX2 .LE. 0.D0) GO TO 3 S2 = (Y(3)-Y(2))/DX2 T = SI + DXI*(SI-S2)/(DXI+DX2) IF (SI .GE. 0.D0) THEN YP(1) = MIN(MAX(0.D0,T), 3.D0*SI) ELSE YP(1) = MAX(MIN(0.D0,T), 3.D0*SI) ENDIF C C YP(I) = (DXIM1*SI+DXI*SIM1)/(DXIM1+DXI) subject to the C constraint that YP(I) has the sign of either SIM1 or C SI, whichever has larger magnitude, and abs(YP(I)) .LE. C 3*min(abs(SIM1),abs(SI)). C DO 1 I = 2,NM1 DXIM1 = DXI DXI = X(I+1) - X(I) IF (DXI .LE. 0.D0) GO TO 3 SIM1 = SI SI = (Y(I+1)-Y(I))/DXI T = (DXIM1*SI+DXI*SIM1)/(DXIM1+DXI) ASIM1 = ABS(SIM1) ASI = ABS(SI) SGN = SIGN(1.D0,SI) IF (ASIM1 .GT. ASI) SGN = SIGN(1.D0,SIM1) IF (SGN .GT. 0.D0) THEN YP(I) = MIN(MAX(0.D0,T),3.D0*MIN(ASIM1,ASI)) ELSE YP(I) = MAX(MIN(0.D0,T),-3.D0*MIN(ASIM1,ASI)) ENDIF 1 CONTINUE C C YP(N) = SNM1 + DXNM1*(SNM1-SNM2)/(DXNM2+DXNM1) subject to C the constraint that YP(N) has the sign of SNM1 and C abs(YP(N)) .LE. 3*abs(SNM1). Note that DXI = DXNM1 and C SI = SNM1. C T = SI + DXI*(SI-SIM1)/(DXIM1+DXI) IF (SI .GE. 0.D0) THEN YP(N) = MIN(MAX(0.D0,T), 3.D0*SI) ELSE YP(N) = MAX(MIN(0.D0,T), 3.D0*SI) ENDIF IER = 0 RETURN C C N is outside its valid range. C 2 IER = 1 RETURN C C X(I+1) .LE. X(I). C 3 IER = I + 1 RETURN END C ---------------------------------------------------------- SUBROUTINE YPC1P (N,X,Y, YP,IER) INTEGER N, IER DOUBLE PRECISION X(N), Y(N), YP(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This subroutine employs a three-point quadratic interpo- C lation method to compute local derivative estimates YP C associated with a set of N data points (X(I),Y(I)). It C is assumed that Y(N) = Y(1), and YP(N) = YP(1) on output. C Thus, a Hermite interpolant H(x) defined by the data C points and derivative estimates is periodic with period C X(N)-X(1). The derivative-estimation formula is the C monotonicity-constrained parabolic fit described in the C reference cited below: H(x) is monotonic in intervals in C which the data is monotonic. The method is invariant C under a linear scaling of the coordinates but is not C additive. C C On input: C C N = Number of data points. N .GE. 3. C C X = Array of length N containing a strictly in- C creasing sequence of abscissae: X(I) < X(I+1) C for I = 1,...,N-1. C C Y = Array of length N containing data values asso- C ciated with the abscissae. Y(N) is set to Y(1) C on output unless IER = 1. C C Input parameters, other than Y(N), are not altered by C this routine. C C On output: C C YP = Array of length N containing estimated deriv- C atives at the abscissae unless IER .NE. 0. C YP is not altered if IER = 1, and is only par- C tially defined if IER > 1. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N < 3. C IER = I if X(I) .LE. X(I-1) for some I in the C range 2,...,N. C C Reference: J. M. Hyman, "Accurate Monotonicity-preserving C Cubic Interpolation", LA-8796-MS, Los C Alamos National Lab, Feb. 1982. C C Modules required by YPC1P: None C C Intrinsic functions called by YPC1P: ABS, MAX, MIN, SIGN C C*********************************************************** C INTEGER I, NM1 DOUBLE PRECISION ASI, ASIM1, DXI, DXIM1, SGN, SI, . SIM1, T C NM1 = N - 1 IF (NM1 .LT. 2) GO TO 2 Y(N) = Y(1) C C Initialize for loop on interior points. C I = 1 DXI = X(2) - X(1) IF (DXI .LE. 0.D0) GO TO 3 SI = (Y(2)-Y(1))/DXI C C YP(I) = (DXIM1*SI+DXI*SIM1)/(DXIM1+DXI) subject to the C constraint that YP(I) has the sign of either SIM1 or C SI, whichever has larger magnitude, and abs(YP(I)) .LE. C 3*min(abs(SIM1),abs(SI)). C DO 1 I = 2,NM1 DXIM1 = DXI DXI = X(I+1) - X(I) IF (DXI .LE. 0.D0) GO TO 3 SIM1 = SI SI = (Y(I+1)-Y(I))/DXI T = (DXIM1*SI+DXI*SIM1)/(DXIM1+DXI) ASIM1 = ABS(SIM1) ASI = ABS(SI) SGN = SIGN(1.D0,SI) IF (ASIM1 .GT. ASI) SGN = SIGN(1.D0,SIM1) IF (SGN .GT. 0.D0) THEN YP(I) = MIN(MAX(0.D0,T),3.D0*MIN(ASIM1,ASI)) ELSE YP(I) = MAX(MIN(0.D0,T),-3.D0*MIN(ASIM1,ASI)) ENDIF 1 CONTINUE C C YP(N) = YP(1), I = 1, and IM1 = N-1. C DXIM1 = DXI DXI = X(2) - X(1) SIM1 = SI SI = (Y(2) - Y(1))/DXI T = (DXIM1*SI + DXI*SIM1)/(DXIM1+DXI) ASIM1 = ABS(SIM1) ASI = ABS(SI) SGN = SIGN(1.D0,SI) IF (ASIM1 .GT. ASI) SGN = SIGN(1.D0,SIM1) IF (SGN .GT. 0.D0) THEN YP(1) = MIN(MAX(0.D0,T), 3.D0*MIN(ASIM1,ASI)) ELSE YP(1) = MAX(MIN(0.D0,T), -3.D0*MIN(ASIM1,ASI)) ENDIF YP(N) = YP(1) C C No error encountered. C IER = 0 RETURN C C N is outside its valid range. C 2 IER = 1 RETURN C C X(I+1) .LE. X(I). C 3 IER = I + 1 RETURN END C ---------------------------------------------------------- SUBROUTINE YPC2 (N,X,Y,SIGMA,ISL1,ISLN,BV1,BVN, . WK, YP,IER) INTEGER N, ISL1, ISLN, IER DOUBLE PRECISION X(N), Y(N), SIGMA(N), BV1, BVN, . WK(N), YP(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This subroutine solves a linear system for a set of C first derivatives YP associated with a Hermite interpola- C tory tension spline H(x). The derivatives are chosen so C that H(x) has two continuous derivatives for all x and H C satisfies user-specified end conditions. C C On input: C C N = Number of data points. N .GE. 2. C C X = Array of length N containing a strictly in- C creasing sequence of abscissae: X(I) < X(I+1) C for I = 1,...,N-1. C C Y = Array of length N containing data values asso- C ciated with the abscissae. H(X(I)) = Y(I) for C I = 1,...,N. C C SIGMA = Array of length N-1 containing tension C factors. SIGMA(I) is associated with inter- C val (X(I),X(I+1)) for I = 1,...,N-1. If C SIGMA(I) = 0, H is the Hermite cubic interp- C olant of the data values and computed deriv- C atives at X(I) and X(I+1), and if all C tension factors are zero, H is the C-2 cubic C spline interpolant which satisfies the end C conditions. C C ISL1 = Option indicator for the condition at X(1): C ISL1 = 0 if YP(1) is to be estimated inter- C nally by a constrained parabolic C fit to the first three points. C This is identical to the method used C by Subroutine YPC1. BV1 is not used C in this case. C ISL1 = 1 if the first derivative of H at X(1) C is specified by BV1. C ISL1 = 2 if the second derivative of H at C X(1) is specified by BV1. C ISL1 = 3 if YP(1) is to be estimated inter- C nally from the derivative of the C tension spline (using SIGMA(1)) C which interpolates the first three C data points and has third derivative C equal to zero at X(1). Refer to C ENDSLP. BV1 is not used in this C case. C C ISLN = Option indicator for the condition at X(N): C ISLN = 0 if YP(N) is to be estimated inter- C nally by a constrained parabolic C fit to the last three data points. C This is identical to the method used C by Subroutine YPC1. BVN is not used C in this case. C ISLN = 1 if the first derivative of H at X(N) C is specified by BVN. C ISLN = 2 if the second derivative of H at C X(N) is specified by BVN. C ISLN = 3 if YP(N) is to be estimated inter- C nally from the derivative of the C tension spline (using SIGMA(N-1)) C which interpolates the last three C data points and has third derivative C equal to zero at X(N). Refer to C ENDSLP. BVN is not used in this C case. C C BV1,BVN = Boundary values or dummy parameters as C defined by ISL1 and ISLN. C C The above parameters are not altered by this routine. C C WK = Array of length at least N-1 to be used as C temporary work space. C C YP = Array of length .GE. N. C C On output: C C YP = Array containing derivatives of H at the C abscissae. YP is not defined if IER .NE. 0. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N, ISL1, or ISLN is outside its C valid range. C IER = I if X(I) .LE. X(I-1) for some I in the C range 2,...,N. C C Modules required by YPC2: ENDSLP, SNHCSH, YPCOEF C C Intrinsic function called by YPC2: ABS C C*********************************************************** C INTEGER I, NM1, NN DOUBLE PRECISION D, D1, D2, DX, R1, R2, S, SD1, SD2, . SIG, YP1, YPN DOUBLE PRECISION ENDSLP C NN = N IF (NN .LT. 2 .OR. ISL1 .LT. 0 .OR. ISL1 .GT. 3 . .OR. ISLN .LT. 0 .OR. ISLN .GT. 3) GO TO 3 NM1 = NN - 1 C C Set YP1 and YPN to the endpoint values. C IF (ISL1 .EQ. 0) THEN IF (NN .GT. 2) YP1 = ENDSLP (X(1),X(2),X(3),Y(1), . Y(2),Y(3),0.D0) ELSEIF (ISL1 .NE. 3) THEN YP1 = BV1 ELSE IF (NN .GT. 2) YP1 = ENDSLP (X(1),X(2),X(3),Y(1), . Y(2),Y(3),SIGMA(1)) ENDIF IF (ISLN .EQ. 0) THEN IF (NN .GT. 2) YPN = ENDSLP (X(NN),X(NM1),X(NN-2), . Y(NN),Y(NM1),Y(NN-2),0.D0) ELSEIF (ISLN .NE. 3) THEN YPN = BVN ELSE IF (NN .GT. 2) YPN = ENDSLP (X(NN),X(NM1),X(NN-2), . Y(NN),Y(NM1),Y(NN-2),SIGMA(NM1)) ENDIF C C Solve the symmetric positive-definite tridiagonal linear C system. The forward elimination step consists of div- C iding each row by its diagonal entry, then introducing a C zero below the diagonal. This requires saving only the C superdiagonal (in WK) and the right hand side (in YP). C I = 1 DX = X(2) - X(1) IF (DX .LE. 0.D0) GO TO 4 S = (Y(2)-Y(1))/DX IF (NN .EQ. 2) THEN IF (ISL1 .EQ. 0 .OR. ISL1 .EQ. 3) YP1 = S IF (ISLN .EQ. 0 .OR. ISLN .EQ. 3) YPN = S ENDIF C C Begin forward elimination. C SIG = ABS(SIGMA(1)) CALL YPCOEF (SIG,DX, D1,SD1) R1 = (SD1+D1)*S WK(1) = 0.D0 YP(1) = YP1 IF (ISL1 .EQ. 2) THEN WK(1) = SD1/D1 YP(1) = (R1-YP1)/D1 ENDIF DO 1 I = 2,NM1 DX = X(I+1) - X(I) IF (DX .LE. 0.D0) GO TO 4 S = (Y(I+1)-Y(I))/DX SIG = ABS(SIGMA(I)) CALL YPCOEF (SIG,DX, D2,SD2) R2 = (SD2+D2)*S D = D1 + D2 - SD1*WK(I-1) WK(I) = SD2/D YP(I) = (R1 + R2 - SD1*YP(I-1))/D D1 = D2 SD1 = SD2 R1 = R2 1 CONTINUE D = D1 - SD1*WK(NM1) YP(NN) = YPN IF (ISLN .EQ. 2) YP(NN) = (R1 + YPN - SD1*YP(NM1))/D C C Back substitution: C DO 2 I = NM1,1,-1 YP(I) = YP(I) - WK(I)*YP(I+1) 2 CONTINUE IER = 0 RETURN C C Invalid integer input parameter. C 3 IER = 1 RETURN C C Abscissae out of order or duplicate points. C 4 IER = I + 1 RETURN END C ---------------------------------------------------------- SUBROUTINE YPC2P (N,X,Y,SIGMA,WK, YP,IER) INTEGER N, IER DOUBLE PRECISION X(N), Y(N), SIGMA(N), WK(*), YP(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This subroutine solves a linear system for a set of C first derivatives YP associated with a Hermite interpola- C tory tension spline H(x). The derivatives are chosen so C that H(x) has two continuous derivatives for all x, and H C satisfies periodic end conditions: first and second der- C ivatives of H at X(1) agree with those at X(N), and thus C the length of a period is X(N) - X(1). It is assumed that C Y(N) = Y(1), and Y(N) is not referenced. C C On input: C C N = Number of data points. N .GE. 3. C C X = Array of length N containing a strictly in- C creasing sequence of abscissae: X(I) < X(I+1) C for I = 1,...,N-1. C C Y = Array of length N containing data values asso- C ciated with the abscissae. H(X(I)) = Y(I) for C I = 1,...,N. C C SIGMA = Array of length N-1 containing tension C factors. SIGMA(I) is associated with inter- C val (X(I),X(I+1)) for I = 1,...,N-1. If C SIGMA(I) = 0, H is the Hermite cubic interp- C olant of the data values and computed deriv- C atives at X(I) and X(I+1), and if all C tension factors are zero, H is the C-2 cubic C spline interpolant which satisfies the end C conditions. C C The above parameters are not altered by this routine. C C WK = Array of length at least 2N-2 to be used as C temporary work space. C C YP = Array of length .GE. N. C C On output: C C YP = Array containing derivatives of H at the C abscissae. YP is not defined if IER .NE. 0. C C IER = Error indicator: C IER = 0 if no errors were encountered. C IER = 1 if N is outside its valid range. C IER = I if X(I) .LE. X(I-1) for some I in the C range 2,...,N. C C Modules required by YPC2P: SNHCSH, YPCOEF C C Intrinsic function called by YPC2P: ABS C C*********************************************************** C INTEGER I, NM1, NM2, NM3, NN, NP1, NPI DOUBLE PRECISION D, D1, D2, DIN, DNM1, DX, R1, R2, . RNM1, S, SD1, SD2, SDNM1, SIG, YPNM1 C NN = N IF (NN .LT. 3) GO TO 4 NM1 = NN - 1 NM2 = NN - 2 NM3 = NN - 3 NP1 = NN + 1 C C The system is order N-1, symmetric, positive-definite, and C tridiagonal except for nonzero elements in the upper C right and lower left corners. The forward elimination C step zeros the subdiagonal and divides each row by its C diagonal entry for the first N-2 rows. The superdiago- C nal is stored in WK(I), the negative of the last column C (fill-in) in WK(N+I), and the right hand side in YP(I) C for I = 1,...,N-2. C I = NM1 DX = X(NN) - X(NM1) IF (DX .LE. 0.D0) GO TO 5 S = (Y(1)-Y(NM1))/DX SIG = ABS(SIGMA(NM1)) CALL YPCOEF (SIG,DX, DNM1,SDNM1) RNM1 = (SDNM1+DNM1)*S I = 1 DX = X(2) - X(1) IF (DX .LE. 0.D0) GO TO 5 S = (Y(2)-Y(1))/DX SIG = ABS(SIGMA(1)) CALL YPCOEF (SIG,DX, D1,SD1) R1 = (SD1+D1)*S D = DNM1 + D1 WK(1) = SD1/D WK(NP1) = -SDNM1/D YP(1) = (RNM1+R1)/D DO 1 I = 2,NM2 DX = X(I+1) - X(I) IF (DX .LE. 0.D0) GO TO 5 S = (Y(I+1)-Y(I))/DX SIG = ABS(SIGMA(I)) CALL YPCOEF (SIG,DX, D2,SD2) R2 = (SD2+D2)*S D = D1 + D2 - SD1*WK(I-1) DIN = 1.D0/D WK(I) = SD2*DIN NPI = NN + I WK(NPI) = -SD1*WK(NPI-1)*DIN YP(I) = (R1 + R2 - SD1*YP(I-1))*DIN SD1 = SD2 D1 = D2 R1 = R2 1 CONTINUE C C The backward elimination step zeros the superdiagonal C (first N-3 elements). WK(I) and YP(I) are overwritten C with the negative of the last column and the new right C hand side, respectively, for I = N-2, N-3, ..., 1. C NPI = NN + NM2 WK(NM2) = WK(NPI) - WK(NM2) DO 2 I = NM3,1,-1 YP(I) = YP(I) - WK(I)*YP(I+1) NPI = NN + I WK(I) = WK(NPI) - WK(I)*WK(I+1) 2 CONTINUE C C Solve the last equation for YP(N-1). C YPNM1 = (R1 + RNM1 - SDNM1*YP(1) - SD1*YP(NM2))/ . (D1 + DNM1 + SDNM1*WK(1) + SD1*WK(NM2)) C C Back substitute for the remainder of the solution C components. C YP(NM1) = YPNM1 DO 3 I = 1,NM2 YP(I) = YP(I) + WK(I)*YPNM1 3 CONTINUE C C YP(N) = YP(1). C YP(N) = YP(1) IER = 0 RETURN C C N is outside its valid range. C 4 IER = 1 RETURN C C Abscissae out of order or duplicate points. C 5 IER = I + 1 RETURN END C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION HPPVAL (T,N,X,Y,YP, . SIGMA, IER) INTEGER N, IER DOUBLE PRECISION T, X(N), Y(N), YP(N), SIGMA(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C This function evaluates the second derivative HPP of a C Hermite interpolatory tension spline H at a point T. C C On input: C C T = Point at which HPP is to be evaluated. Extrap- C olation is performed if T < X(1) or T > X(N). C C N = Number of data points. N .GE. 2. C C X = Array of length N containing the abscissae. C These must be in strictly increasing order: C X(I) < X(I+1) for I = 1,...,N-1. C C Y = Array of length N containing data values. C H(X(I)) = Y(I) for I = 1,...,N. C C YP = Array of length N containing first deriva- C tives. HP(X(I)) = YP(I) for I = 1,...,N, where C HP denotes the derivative of H. C C SIGMA = Array of length N-1 containing tension fac- C tors whose absolute values determine the C balance between cubic and linear in each C interval. SIGMA(I) is associated with int- C erval (I,I+1) for I = 1,...,N-1. C C Input parameters are not altered by this function. C C On output: C C IER = Error indicator: C IER = 0 if no errors were encountered and C X(1) .LE. T .LE. X(N). C IER = 1 if no errors were encountered and C extrapolation was necessary. C IER = -1 if N < 2. C IER = -2 if the abscissae are not in strictly C increasing order. (This error will C not necessarily be detected.) C C HPPVAL = Second derivative value HPP(T), or zero if C IER < 0. C C Modules required by HPPVAL: INTRVL, SNHCSH C C Intrinsic functions called by HPPVAL: ABS, EXP C C*********************************************************** C INTEGER I, IP1 DOUBLE PRECISION B1, B2, CM, CM2, CMM, COSH2, D1, D2, . DUMMY, DX, E, E1, E2, EMS, S, SB1, . SB2, SBIG, SIG, SINH2, SM, SM2, TM INTEGER INTRVL C DATA SBIG/85.D0/ IF (N .LT. 2) GO TO 1 C C Find the index of the left end of an interval containing C T. If T < X(1) or T > X(N), extrapolation is performed C using the leftmost or rightmost interval. C IF (T .LT. X(1)) THEN I = 1 IER = 1 ELSEIF (T .GT. X(N)) THEN I = N-1 IER = 1 ELSE I = INTRVL (T,N,X) IER = 0 ENDIF IP1 = I + 1 C C Compute interval width DX, local coordinates B1 and B2, C and second differences D1 and D2. C DX = X(IP1) - X(I) IF (DX .LE. 0.D0) GO TO 2 B1 = (X(IP1) - T)/DX B2 = 1.D0 - B1 S = (Y(IP1)-Y(I))/DX D1 = S - YP(I) D2 = YP(IP1) - S SIG = ABS(SIGMA(I)) IF (SIG .LT. 1.D-9) THEN C C SIG = 0: H is the Hermite cubic interpolant. C HPPVAL = (D1 + D2 + 3.D0*(B2-B1)*(D2-D1))/DX ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5: use approximations designed to avoid C cancellation error in the hyperbolic functions. C SB2 = SIG*B2 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB2, SM2,CM2,DUMMY) SINH2 = SM2 + SB2 COSH2 = CM2 + 1.D0 E = SIG*SM - CMM - CMM HPPVAL = SIG*((CM*SINH2-SM*COSH2)*(D1+D2) + . SIG*(CM*COSH2-(SM+SIG)*SINH2)*D1)/(DX*E) ELSE C C SIG > .5: use negative exponentials in order to avoid C overflow. Note that EMS = EXP(-SIG). In the case of C extrapolation (negative B1 or B2), H is approximated by C a linear function if -SIG*B1 or -SIG*B2 is large. C SB1 = SIG*B1 SB2 = SIG - SB1 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN HPPVAL = 0.D0 ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1.D0 - EMS E = TM*(SIG*(1.D0+EMS) - TM - TM) HPPVAL = SIG*(SIG*((E1*EMS+E2)*D1+(E1+E2*EMS)*D2)- . TM*(E1+E2)*(D1+D2))/(DX*E) ENDIF ENDIF RETURN C C N is outside its valid range. C 1 HPPVAL = 0.D0 IER = -1 RETURN C C X(I) .GE. X(I+1). C 2 HPPVAL = 0.D0 IER = -2 RETURN END C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION HPVAL (T,N,X,Y,YP, . SIGMA, IER) INTEGER N, IER DOUBLE PRECISION T, X(N), Y(N), YP(N), SIGMA(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C This function evaluates the first derivative HP of a C Hermite interpolatory tension spline H at a point T. C C On input: C C T = Point at which HP is to be evaluated. Extrapo- C lation is performed if T < X(1) or T > X(N). C C N = Number of data points. N .GE. 2. C C X = Array of length N containing the abscissae. C These must be in strictly increasing order: C X(I) < X(I+1) for I = 1,...,N-1. C C Y = Array of length N containing data values. C H(X(I)) = Y(I) for I = 1,...,N. C C YP = Array of length N containing first deriva- C tives. HP(X(I)) = YP(I) for I = 1,...,N. C C SIGMA = Array of length N-1 containing tension fac- C tors whose absolute values determine the C balance between cubic and linear in each C interval. SIGMA(I) is associated with int- C erval (I,I+1) for I = 1,...,N-1. C C Input parameters are not altered by this function. C C On output: C C IER = Error indicator: C IER = 0 if no errors were encountered and C X(1) .LE. T .LE. X(N). C IER = 1 if no errors were encountered and C extrapolation was necessary. C IER = -1 if N < 2. C IER = -2 if the abscissae are not in strictly C increasing order. (This error will C not necessarily be detected.) C C HPVAL = Derivative value HP(T), or zero if IER < 0. C C Modules required by HPVAL: INTRVL, SNHCSH C C Intrinsic functions called by HPVAL: ABS, EXP C C*********************************************************** C INTEGER I, IP1 DOUBLE PRECISION B1, B2, CM, CM2, CMM, D1, D2, DUMMY, . DX, E, E1, E2, EMS, S, S1, SB1, SB2, . SBIG, SIG, SINH2, SM, SM2, TM INTEGER INTRVL C DATA SBIG/85.D0/ IF (N .LT. 2) GO TO 1 C C Find the index of the left end of an interval containing C T. If T < X(1) or T > X(N), extrapolation is performed C using the leftmost or rightmost interval. C IF (T .LT. X(1)) THEN I = 1 IER = 1 ELSEIF (T .GT. X(N)) THEN I = N-1 IER = 1 ELSE I = INTRVL (T,N,X) IER = 0 ENDIF IP1 = I + 1 C C Compute interval width DX, local coordinates B1 and B2, C and second differences D1 and D2. C DX = X(IP1) - X(I) IF (DX .LE. 0.D0) GO TO 2 B1 = (X(IP1) - T)/DX B2 = 1.D0 - B1 S1 = YP(I) S = (Y(IP1)-Y(I))/DX D1 = S - S1 D2 = YP(IP1) - S SIG = ABS(SIGMA(I)) IF (SIG .LT. 1.D-9) THEN C C SIG = 0: H is the Hermite cubic interpolant. C HPVAL = S1 + B2*(D1 + D2 - 3.D0*B1*(D2-D1)) ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5: use approximations designed to avoid C cancellation error in the hyperbolic functions. C SB2 = SIG*B2 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB2, SM2,CM2,DUMMY) SINH2 = SM2 + SB2 E = SIG*SM - CMM - CMM HPVAL = S1 + ((CM*CM2-SM*SINH2)*(D1+D2) + . SIG*(CM*SINH2-(SM+SIG)*CM2)*D1)/E ELSE C C SIG > .5: use negative exponentials in order to avoid C overflow. Note that EMS = EXP(-SIG). In the case of C extrapolation (negative B1 or B2), H is approximated by C a linear function if -SIG*B1 or -SIG*B2 is large. C SB1 = SIG*B1 SB2 = SIG - SB1 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN HPVAL = S ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1.D0 - EMS E = TM*(SIG*(1.D0+EMS) - TM - TM) HPVAL = S + (TM*((E2-E1)*(D1+D2) + TM*(D1-D2)) + . SIG*((E1*EMS-E2)*D1 + (E1-E2*EMS)*D2))/E ENDIF ENDIF RETURN C C N is outside its valid range. C 1 HPVAL = 0.D0 IER = -1 RETURN C C X(I) .GE. X(I+1). C 2 HPVAL = 0.D0 IER = -2 RETURN END C ---------------------------------------------------------- DOUBLE PRECISION FUNCTION HVAL (T,N,X,Y,YP,SIGMA, IER) INTEGER N, IER DOUBLE PRECISION T, X(N), Y(N), YP(N), SIGMA(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 11/17/96 C C This function evaluates a Hermite interpolatory tension C spline H at a point T. Note that a large value of SIGMA C may cause underflow. The result is assumed to be zero. C C Given arrays X, Y, YP, and SIGMA of length NN, if T is C known to lie in the interval (X(I),X(J)) for some I < J, C a gain in efficiency can be achieved by calling this C function with N = J+1-I (rather than NN) and the I-th C components of the arrays (rather than the first) as par- C ameters. C C On input: C C T = Point at which H is to be evaluated. Extrapo- C lation is performed if T < X(1) or T > X(N). C (TE(I),N,T,X,XP,SIGMA, IERR) C N = Number of data points. N .GE. 2. C C X = Array of length N containing the abscissae. C These must be in strictly increasing order: C X(I) < X(I+1) for I = 1,...,N-1. C C Y = Array of length N containing data values. C H(X(I)) = Y(I) for I = 1,...,N. C C YP = Array of length N containing first deriva- C tives. HP(X(I)) = YP(I) for I = 1,...,N, where C HP denotes the derivative of H. C C SIGMA = Array of length N-1 containing tension fac- C tors whose absolute values determine the C balance between cubic and linear in each C interval. SIGMA(I) is associated with int- C erval (I,I+1) for I = 1,...,N-1. C C Input parameters are not altered by this function. C C On output: C C IER = Error indicator: C IER = 0 if no errors were encountered and C X(1) .LE. T .LE. X(N). C IER = 1 if no errors were encountered and C extrapolation was necessary. C IER = -1 if N < 2. C IER = -2 if the abscissae are not in strictly C increasing order. (This error will C not necessarily be detected.) C C HVAL = Function value H(T), or zero if IER < 0. C C Modules required by HVAL: INTRVL, SNHCSH C C Intrinsic functions called by HVAL: ABS, EXP C C*********************************************************** C INTEGER I, IP1 DOUBLE PRECISION B1, B2, CM, CM2, CMM, D1, D2, DUMMY, . DX, E, E1, E2, EMS, S, S1, SB1, SB2, . SBIG, SIG, SM, SM2, TM, TP, TS, U, Y1 INTEGER INTRVL C DATA SBIG/85.D0/ IF (N .LT. 2) GO TO 1 C C Find the index of the left end of an interval containing C T. If T < X(1) or T > X(N), extrapolation is performed C using the leftmost or rightmost interval. C IF (T .LT. X(1)) THEN I = 1 IER = 1 ELSEIF (T .GT. X(N)) THEN I = N-1 IER = 1 ELSE I = INTRVL (T,N,X) IER = 0 ENDIF IP1 = I + 1 C C Compute interval width DX, local coordinates B1 and B2, C and second differences D1 and D2. C DX = X(IP1) - X(I) IF (DX .LE. 0.D0) GO TO 2 U = T - X(I) B2 = U/DX B1 = 1.D0 - B2 Y1 = Y(I) S1 = YP(I) S = (Y(IP1)-Y1)/DX D1 = S - S1 D2 = YP(IP1) - S SIG = ABS(SIGMA(I)) IF (SIG .LT. 1.D-9) THEN C C SIG = 0: H is the Hermite cubic interpolant. C HVAL = Y1 + U*(S1 + B2*(D1 + B1*(D1-D2))) ELSEIF (SIG .LE. .5D0) THEN C C 0 .LT. SIG .LE. .5: use approximations designed to avoid C cancellation error in the hyperbolic functions. C SB2 = SIG*B2 CALL SNHCSH (SIG, SM,CM,CMM) CALL SNHCSH (SB2, SM2,CM2,DUMMY) E = SIG*SM - CMM - CMM HVAL = Y1 + S1*U + DX*((CM*SM2-SM*CM2)*(D1+D2) + . SIG*(CM*CM2-(SM+SIG)*SM2)*D1) . /(SIG*E) ELSE C C SIG > .5: use negative exponentials in order to avoid C overflow. Note that EMS = EXP(-SIG). In the case of C extrapolation (negative B1 or B2), H is approximated by C a linear function if -SIG*B1 or -SIG*B2 is large. C SB1 = SIG*B1 SB2 = SIG - SB1 IF (-SB1 .GT. SBIG .OR. -SB2 .GT. SBIG) THEN HVAL = Y1 + S*U ELSE E1 = EXP(-SB1) E2 = EXP(-SB2) EMS = E1*E2 TM = 1.D0 - EMS TS = TM*TM TP = 1.D0 + EMS E = TM*(SIG*TP - TM - TM) HVAL = Y1 + S*U + DX*(TM*(TP-E1-E2)*(D1+D2) + SIG* . ((E2+EMS*(E1-2.D0)-B1*TS)*D1+ . (E1+EMS*(E2-2.D0)-B2*TS)*D2))/ . (SIG*E) ENDIF ENDIF RETURN C C N is outside its valid range. C 1 HVAL = 0.D0 IER = -1 RETURN C C X(I) .GE. X(I+1). C 2 HVAL = 0.D0 IER = -2 RETURN END C ---------------------------------------------------------- INTEGER FUNCTION INTRVL (T,N,X) INTEGER N DOUBLE PRECISION T, X(N) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This function returns the index of the left end of an C interval (defined by an increasing sequence X) which C contains the value T. The method consists of first test- C ing the interval returned by a previous call, if any, and C then using a binary search if necessary. C C On input: C C T = Point to be located. C C N = Length of X. N .GE. 2. C C X = Array of length N assumed (without a test) to C contain a strictly increasing sequence of C values. C C Input parameters are not altered by this function. C C On output: C C INTRVL = Index I defined as follows: C C I = 1 if T .LT. X(2) or N .LE. 2, C I = N-1 if T .GE. X(N-1), and C X(I) .LE. T .LT. X(I+1) otherwise. C C Modules required by INTRVL: None C C*********************************************************** C INTEGER IH, IL, K DOUBLE PRECISION TT C SAVE IL DATA IL/1/ TT = T IF (IL .GE. 1 .AND. IL .LT. N) THEN IF (X(IL) .LE. TT .AND. TT .LT. X(IL+1)) GO TO 2 ENDIF C C Initialize low and high indexes. C IL = 1 IH = N C C Binary search: C 1 IF (IH .LE. IL+1) GO TO 2 K = (IL+IH)/2 IF (TT .LT. X(K)) THEN IH = K ELSE IL = K ENDIF GO TO 1 C C X(IL) .LE. T .LT. X(IL+1) or (T .LT. X(1) and IL=1) C or (T .GE. X(N) and IL=N-1) C 2 INTRVL = IL RETURN END C ---------------------------------------------------------- SUBROUTINE TSVAL1 (N,X,Y,YP,SIGMA,IFLAG,NE,TE, V, . IER) INTEGER N, IFLAG, NE, IER DOUBLE PRECISION X(N), Y(N), YP(N), SIGMA(N), TE(NE), . V(NE) C C*********************************************************** C C From TSPACK C Robert J. Renka C Dept. of Computer Science C Univ. of North Texas C renka@cs.unt.edu C 06/10/92 C C This subroutine evaluates a Hermite interpolatory ten- C sion spline H or its first or second derivative at a set C of points TE. C C Note that a large tension factor in SIGMA may cause C underflow. The result is assumed to be zero. If not the C default, this may be specified by either a compiler option C or operating system option. C C On input: C C N = Number of data points. N .GE. 2. C C X = Array of length N containing the abscissae. C These must be in strictly increasing order: C X(I) < X(I+1) for I = 1,...,N-1. C C Y = Array of length N containing data values or C function values returned by Subroutine SMCRV. C Y(I) = H(X(I)) for I = 1,...,N. C C YP = Array of length N containing first deriva- C tives. YP(I) = HP(X(I)) for I = 1,...,N, where C HP denotes the derivative of H. C C SIGMA = Array of length N-1 containing tension fac- C tors whose absolute values determine the C balance between cubic and linear in each C interval. SIGMA(I) is associated with int- C erval (I,I+1) for I = 1,...,N-1. C C IFLAG = Output option indicator: C IFLAG = 0 if values of H are to be computed. C IFLAG = 1 if first derivative values are to C be computed. C IFLAG = 2 if second derivative values are to C be computed. C C NE = Number of evaluation points. NE > 0. C C TE = Array of length NE containing the evaluation C points. The sequence should be strictly in- C creasing for maximum efficiency. Extrapolation C is performed if a point is not in the interval C [X(1),X(N)]. C C The above parameters are not altered by this routine. C C V = Array of length at least NE. C C On output: C C V = Array of function, first derivative, or second C derivative values at the evaluation points un- C less IER < 0. If IER = -1, V is not altered. C If IER = -2, V may be only partially defined. C C IER = Error indicator: C IER = 0 if no errors were encountered and C no extrapolation occurred. C IER > 0 if no errors were encountered but C extrapolation was required at IER C points. C IER = -1 if N < 2, IFLAG < 0, IFLAG > 2, or C NE < 1. C IER = -2 if the abscissae are not in strictly C increasing order. (This error will C not necessarily be detected.) C C Modules required by TSVAL1: HPPVAL, HPVAL, HVAL, INTRVL, C SNHCSH C C*********************************************************** C INTEGER I, IERR, IFLG, NVAL, NX DOUBLE PRECISION HPPVAL, HPVAL, HVAL C IFLG = IFLAG NVAL = NE C C Test for invalid input. C IF (N .LT. 2 .OR. IFLG .LT. 0 .OR. IFLG .GT. 2 . .OR. NVAL .LT. 1) GO TO 2 C C Initialize the number of extrapolation points NX and C loop on evaluation points. C NX = 0 DO 1 I = 1,NVAL IF (IFLG .EQ. 0) THEN V(I) = HVAL (TE(I),N,X,Y,YP,SIGMA, IERR) ELSEIF (IFLG .EQ. 1) THEN V(I) = HPVAL (TE(I),N,X,Y,YP,SIGMA, IERR) ELSE V(I) = HPPVAL (TE(I),N,X,Y,YP,SIGMA, IERR) ENDIF IF (IERR .LT. 0) GO TO 3 NX = NX + IERR 1 CONTINUE C C No errors encountered. C IER = NX RETURN C C N, IFLAG, or NE is outside its valid range. C 2 IER = -1 RETURN C C X is not strictly increasing. C 3 IER = -2 RETURN END C*********************************************************** C Minpack Copyright Notice (1999) University of Chicago. All rights reserved C Redistribution and use in source and binary forms, with or C without modification, are permitted provided that the C following conditions are met: C 1. Redistributions of source code must retain the above C copyright notice, this list of conditions and the following C disclaimer. C 2. Redistributions in binary form must reproduce the above C copyright notice, this list of conditions and the following C disclaimer in the documentation and/or other materials C provided with the distribution. C 3. The end-user documentation included with the C redistribution, if any, must include the following C acknowledgment: C "This product includes software developed by the C University of Chicago, as Operator of Argonne National C Laboratory. C Alternately, this acknowledgment may appear in the software C itself, if and wherever such third-party acknowledgments C normally appear. C 4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" C WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE C UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND C THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR C IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES C OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE C OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY C OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR C USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF C THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) C DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION C UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL C BE CORRECTED. C 5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT C HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF C ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, C INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF C ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF C PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER C SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT C (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, C EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE C POSSIBILITY OF SUCH LOSS OR DAMAGES. C*********************************************************** SUBROUTINE FCN(m,n,x,FVEC,IFLAG,xPSF,yPSF) INTEGER i,m,n,IFLAG DOUBLE PRECISION FVEC(m),x(m) DOUBLE PRECISION Mof(m),xPSF(m),yPSF(m) DO i = 1,m Mof(i) = (1.0D0 + (xPSF(i)/x(1))**2.D0)**(-1.D0*x(2)) FVEC(i) = yPSF(i) - Mof(i) ENDDO END C*********************************************************** double precision function enorm(n,x) integer n double precision x(n) c ********** integer i double precision agiant,floatn,one,rdwarf,rgiant,s1,s2,s3,xabs, * x1max,x3max,zero data one,zero,rdwarf,rgiant /1.0d0,0.0d0,3.834d-20,1.304d19/ s1 = zero s2 = zero s3 = zero x1max = zero x3max = zero floatn = n agiant = rgiant/floatn do 90 i = 1, n xabs = dabs(x(i)) if (xabs .gt. rdwarf .and. xabs .lt. agiant) go to 70 if (xabs .le. rdwarf) go to 30 c c sum for large components. c if (xabs .le. x1max) go to 10 s1 = one + s1*(x1max/xabs)**2 x1max = xabs go to 20 10 continue s1 = s1 + (xabs/x1max)**2 20 continue go to 60 30 continue c c sum for small components. c if (xabs .le. x3max) go to 40 s3 = one + s3*(x3max/xabs)**2 x3max = xabs go to 50 40 continue if (xabs .ne. zero) s3 = s3 + (xabs/x3max)**2 50 continue 60 continue go to 80 70 continue c c sum for intermediate components. c s2 = s2 + xabs**2 80 continue 90 continue c c calculation of norm. c if (s1 .eq. zero) go to 100 enorm = x1max*dsqrt(s1+(s2/x1max)/x1max) go to 130 100 continue if (s2 .eq. zero) go to 110 if (s2 .ge. x3max) * enorm = dsqrt(s2*(one+(x3max/s2)*(x3max*s3))) if (s2 .lt. x3max) * enorm = dsqrt(x3max*((s2/x3max)+(x3max*s3))) go to 120 110 continue enorm = x3max*dsqrt(s3) 120 continue 130 continue return c c last card of function enorm. c end C*********************************************************** double precision function dpmpar(i) integer i c ********** integer mcheps(4) integer minmag(4) integer maxmag(4) double precision dmach(3) equivalence (dmach(1),mcheps(1)) equivalence (dmach(2),minmag(1)) equivalence (dmach(3),maxmag(1)) c c c Machine constants for IEEE machines. c data dmach(1) /2.22044604926d-16/ data dmach(2) /2.22507385852d-308/ data dmach(3) /1.79769313485d+308/ c dpmpar = dmach(i) return c end C*********************************************************** subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa &,xPSF,yPSF) integer m,n,ldfjac,iflag double precision epsfcn,xPSF(m),yPSF(m) double precision x(n),fvec(m),fjac(ldfjac,n),wa(m) c ********** integer i,j double precision eps,epsmch,h,temp,zero double precision dpmpar data zero /0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c eps = dsqrt(dmax1(epsfcn,epsmch)) do 20 j = 1, n temp = x(j) h = eps*dabs(temp) if (h .eq. zero) h = eps x(j) = temp + h call fcn(m,n,x,wa,iflag,xPSF,yPSF) if (iflag .lt. 0) go to 30 x(j) = temp do 10 i = 1, m fjac(i,j) = (wa(i) - fvec(i))/h 10 continue 20 continue 30 continue return c c end C*********************************************************** subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn, * diag,mode,factor,nprint,info,nfev,fjac,ldfjac, * ipvt,qtf,wa1,wa2,wa3,wa4,xPSF,yPSF) integer m,n,maxfev,mode,nprint,info,nfev,ldfjac integer ipvt(n) double precision ftol,xtol,gtol,epsfcn,factor double precision x(n),fvec(m),diag(n),fjac(ldfjac,n),qtf(n), * wa1(n),wa2(n),wa3(n),wa4(m),xPSF(m),yPSF(m) external fcn c ********** c c subroutine lmdif c c the purpose of lmdif is to minimize the sum of the squares of c m nonlinear functions in n variables by a modification of c the levenberg-marquardt algorithm. the user must provide a c subroutine which calculates the functions. the jacobian is c then calculated by a forward-difference approximation. c c the subroutine statement is c c subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn, c diag,mode,factor,nprint,info,nfev,fjac, c ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4) c c where c c fcn is the name of the user-supplied subroutine which c calculates the functions. fcn must be declared c in an external statement in the user calling c program, and should be written as follows. c c subroutine fcn(m,n,x,fvec,iflag) c integer m,n,iflag c double precision x(n),fvec(m) c ---------- c calculate the functions at x and c return this vector in fvec. c ---------- c return c end c c the value of iflag should not be changed by fcn unless c the user wants to terminate execution of lmdif. c in this case set iflag to a negative integer. c c m is a positive integer input variable set to the number c of functions. c c n is a positive integer input variable set to the number c of variables. n must not exceed m. c c x is an array of length n. on input x must contain c an initial estimate of the solution vector. on output x c contains the final estimate of the solution vector. c c fvec is an output array of length m which contains c the functions evaluated at the output x. c c ftol is a nonnegative input variable. termination c occurs when both the actual and predicted relative c reductions in the sum of squares are at most ftol. c therefore, ftol measures the relative error desired c in the sum of squares. c c xtol is a nonnegative input variable. termination c occurs when the relative error between two consecutive c iterates is at most xtol. therefore, xtol measures the c relative error desired in the approximate solution. c c gtol is a nonnegative input variable. termination c occurs when the cosine of the angle between fvec and c any column of the jacobian is at most gtol in absolute c value. therefore, gtol measures the orthogonality c desired between the function vector and the columns c of the jacobian. c c maxfev is a positive integer input variable. termination c occurs when the number of calls to fcn is at least c maxfev by the end of an iteration. c c epsfcn is an input variable used in determining a suitable c step length for the forward-difference approximation. this c approximation assumes that the relative errors in the c functions are of the order of epsfcn. if epsfcn is less c than the machine precision, it is assumed that the relative c errors in the functions are of the order of the machine c precision. c c diag is an array of length n. if mode = 1 (see c below), diag is internally set. if mode = 2, diag c must contain positive entries that serve as c multiplicative scale factors for the variables. c c mode is an integer input variable. if mode = 1, the c variables will be scaled internally. if mode = 2, c the scaling is specified by the input diag. other c values of mode are equivalent to mode = 1. c c factor is a positive input variable used in determining the c initial step bound. this bound is set to the product of c factor and the euclidean norm of diag*x if nonzero, or else c to factor itself. in most cases factor should lie in the c interval (.1,100.). 100. is a generally recommended value. c c nprint is an integer input variable that enables controlled c write(*,*)ing of iterates if it is positive. in this case, c fcn is called with iflag = 0 at the beginning of the first c iteration and every nprint iterations thereafter and c immediately prior to return, with x and fvec available c for write(*,*)ing. if nprint is not positive, no special calls c of fcn with iflag = 0 are made. c c info is an integer output variable. if the user has c terminated execution, info is set to the (negative) c value of iflag. see description of fcn. otherwise, c info is set as follows. c c info = 0 improper input parameters. c c info = 1 both actual and predicted relative reductions c in the sum of squares are at most ftol. c c info = 2 relative error between two consecutive iterates c is at most xtol. c c info = 3 conditions for info = 1 and info = 2 both hold. c c info = 4 the cosine of the angle between fvec and any c column of the jacobian is at most gtol in c absolute value. c c info = 5 number of calls to fcn has reached or c exceeded maxfev. c c info = 6 ftol is too small. no further reduction in c the sum of squares is possible. c c info = 7 xtol is too small. no further improvement in c the approximate solution x is possible. c c info = 8 gtol is too small. fvec is orthogonal to the c columns of the jacobian to machine precision. c c nfev is an integer output variable set to the number of c calls to fcn. c c fjac is an output m by n array. the upper n by n submatrix c of fjac contains an upper triangular matrix r with c diagonal elements of nonincreasing magnitude such that c c t t t c p *(jac *jac)*p = r *r, c c where p is a permutation matrix and jac is the final c calculated jacobian. column j of p is column ipvt(j) c (see below) of the identity matrix. the lower trapezoidal c part of fjac contains information generated during c the computation of r. c c ldfjac is a positive integer input variable not less than m c which specifies the leading dimension of the array fjac. c c ipvt is an integer output array of length n. ipvt c defines a permutation matrix p such that jac*p = q*r, c where jac is the final calculated jacobian, q is c orthogonal (not stored), and r is upper triangular c with diagonal elements of nonincreasing magnitude. c column j of p is column ipvt(j) of the identity matrix. c c qtf is an output array of length n which contains c the first n elements of the vector (q transpose)*fvec. c c wa1, wa2, and wa3 are work arrays of length n. c c wa4 is a work array of length m. c c subprograms called c c user-supplied ...... fcn c c minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac c c fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod c c argonne national laboratory. minpack project. march 1980. c burton s. garbow, kenneth e. hillstrom, jorge j. more c c ********** integer i,iflag,iter,j,l double precision actred,delta,dirder,epsmch,fnorm,fnorm1,gnorm, * one,par,pnorm,prered,p1,p5,p25,p75,p0001,ratio, * sum,temp,temp1,temp2,xnorm,zero double precision dpmpar,enorm data one,p1,p5,p25,p75,p0001,zero * /1.0d0,1.0d-1,5.0d-1,2.5d-1,7.5d-1,1.0d-4,0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c info = 0 iflag = 0 nfev = 0 c c check the input parameters for errors. c if (n .le. 0 .or. m .lt. n .or. ldfjac .lt. m * .or. ftol .lt. zero .or. xtol .lt. zero .or. gtol .lt. zero * .or. maxfev .le. 0 .or. factor .le. zero) go to 300 if (mode .ne. 2) go to 20 do 10 j = 1, n if (diag(j) .le. zero) go to 300 10 continue 20 continue c c evaluate the function at the starting point c and calculate its norm. c iflag = 1 call fcn(m,n,x,fvec,iflag,xPSF,yPSF) nfev = 1 if (iflag .lt. 0) go to 300 fnorm = enorm(m,fvec) c c initialize levenberg-marquardt parameter and iteration counter. c par = zero iter = 1 c c beginning of the outer loop. c 30 continue c c calculate the jacobian matrix. c iflag = 2 call fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa4, &xPSF,yPSF) nfev = nfev + n if (iflag .lt. 0) go to 300 c c if requested, call fcn to enable write(*,*)ing of iterates. c if (nprint .le. 0) go to 40 iflag = 0 if (mod(iter-1,nprint) .eq. 0) call fcn(m,n,x,fvec,iflag, &xPSF,yPSF) if (iflag .lt. 0) go to 300 40 continue c c compute the qr factorization of the jacobian. c call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3) c c on the first iteration and if mode is 1, scale according c to the norms of the columns of the initial jacobian. c if (iter .ne. 1) go to 80 if (mode .eq. 2) go to 60 do 50 j = 1, n diag(j) = wa2(j) if (wa2(j) .eq. zero) diag(j) = one 50 continue 60 continue c c on the first iteration, calculate the norm of the scaled x c and initialize the step bound delta. c do 70 j = 1, n wa3(j) = diag(j)*x(j) 70 continue xnorm = enorm(n,wa3) delta = factor*xnorm if (delta .eq. zero) delta = factor 80 continue c c form (q transpose)*fvec and store the first n components in c qtf. c do 90 i = 1, m wa4(i) = fvec(i) 90 continue do 130 j = 1, n if (fjac(j,j) .eq. zero) go to 120 sum = zero do 100 i = j, m sum = sum + fjac(i,j)*wa4(i) 100 continue temp = -sum/fjac(j,j) do 110 i = j, m wa4(i) = wa4(i) + fjac(i,j)*temp 110 continue 120 continue fjac(j,j) = wa1(j) qtf(j) = wa4(j) 130 continue c c compute the norm of the scaled gradient. c gnorm = zero if (fnorm .eq. zero) go to 170 do 160 j = 1, n l = ipvt(j) if (wa2(l) .eq. zero) go to 150 sum = zero do 140 i = 1, j sum = sum + fjac(i,j)*(qtf(i)/fnorm) 140 continue gnorm = dmax1(gnorm,dabs(sum/wa2(l))) 150 continue 160 continue 170 continue c c test for convergence of the gradient norm. c if (gnorm .le. gtol) info = 4 if (info .ne. 0) go to 300 c c rescale if necessary. c if (mode .eq. 2) go to 190 do 180 j = 1, n diag(j) = dmax1(diag(j),wa2(j)) 180 continue 190 continue c c beginning of the inner loop. c 200 continue c c determine the levenberg-marquardt parameter. c call lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,par,wa1,wa2, * wa3,wa4) c c store the direction p and x + p. calculate the norm of p. c do 210 j = 1, n wa1(j) = -wa1(j) wa2(j) = x(j) + wa1(j) wa3(j) = diag(j)*wa1(j) 210 continue pnorm = enorm(n,wa3) c c on the first iteration, adjust the initial step bound. c if (iter .eq. 1) delta = dmin1(delta,pnorm) c c evaluate the function at x + p and calculate its norm. c iflag = 1 call fcn(m,n,wa2,wa4,iflag,xPSF,yPSF) nfev = nfev + 1 if (iflag .lt. 0) go to 300 fnorm1 = enorm(m,wa4) c c compute the scaled actual reduction. c actred = -one if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2 c c compute the scaled predicted reduction and c the scaled directional derivative. c do 230 j = 1, n wa3(j) = zero l = ipvt(j) temp = wa1(l) do 220 i = 1, j wa3(i) = wa3(i) + fjac(i,j)*temp 220 continue 230 continue temp1 = enorm(n,wa3)/fnorm temp2 = (dsqrt(par)*pnorm)/fnorm prered = temp1**2 + temp2**2/p5 dirder = -(temp1**2 + temp2**2) c c compute the ratio of the actual to the predicted c reduction. c ratio = zero if (prered .ne. zero) ratio = actred/prered c c update the step bound. c if (ratio .gt. p25) go to 240 if (actred .ge. zero) temp = p5 if (actred .lt. zero) * temp = p5*dirder/(dirder + p5*actred) if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1 delta = temp*dmin1(delta,pnorm/p1) par = par/temp go to 260 240 continue if (par .ne. zero .and. ratio .lt. p75) go to 250 delta = pnorm/p5 par = p5*par 250 continue 260 continue c c test for successful iteration. c if (ratio .lt. p0001) go to 290 c c successful iteration. update x, fvec, and their norms. c do 270 j = 1, n x(j) = wa2(j) wa2(j) = diag(j)*x(j) 270 continue do 280 i = 1, m fvec(i) = wa4(i) 280 continue xnorm = enorm(n,wa2) fnorm = fnorm1 iter = iter + 1 290 continue c c tests for convergence. c if (dabs(actred) .le. ftol .and. prered .le. ftol * .and. p5*ratio .le. one) info = 1 if (delta .le. xtol*xnorm) info = 2 if (dabs(actred) .le. ftol .and. prered .le. ftol * .and. p5*ratio .le. one .and. info .eq. 2) info = 3 if (info .ne. 0) go to 300 c c tests for termination and stringent tolerances. c if (nfev .ge. maxfev) info = 5 if (dabs(actred) .le. epsmch .and. prered .le. epsmch * .and. p5*ratio .le. one) info = 6 if (delta .le. epsmch*xnorm) info = 7 if (gnorm .le. epsmch) info = 8 if (info .ne. 0) go to 300 c c end of the inner loop. repeat if iteration unsuccessful. c if (ratio .lt. p0001) go to 200 c c end of the outer loop. c go to 30 300 continue c c termination, either normal or user imposed. c if (iflag .lt. 0) info = iflag iflag = 0 if (nprint .gt. 0) call fcn(m,n,x,fvec,iflag,xPSF,yPSF) return c c end C*********************************************************** subroutine lmdif1(fcn,m,n,x,fvec,tol,info,iwa,wa,lwa,xPSF,yPSF) integer m,n,info,lwa integer iwa(n) double precision tol double precision x(n),fvec(m),wa(lwa),xPSF(m),yPSF(m) external fcn c ********** integer maxfev,mode,mp5n,nfev,nprint double precision epsfcn,factor,ftol,gtol,xtol,zero data factor,zero /1.0d2,0.0d0/ info = 0 c c check the input parameters for errors. c if (n .le. 0 .or. m .lt. n .or. tol .lt. zero * .or. lwa .lt. m*n + 5*n + m) go to 10 c c call lmdif. c maxfev = 200*(n + 1) ftol = tol xtol = tol gtol = zero epsfcn = zero mode = 1 nprint = 0 mp5n = m + 5*n call lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,wa(1), * mode,factor,nprint,info,nfev,wa(mp5n+1),m,iwa, * wa(n+1),wa(2*n+1),wa(3*n+1),wa(4*n+1),wa(5*n+1), * xPSF,yPSF) if (info .eq. 8) info = 4 10 continue return c c end C*********************************************************** subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,wa1, * wa2) integer n,ldr integer ipvt(n) double precision delta,par double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa1(n), * wa2(n) c ********** integer i,iter,j,jm1,jp1,k,l,nsing double precision dxnorm,dwarf,fp,gnorm,parc,parl,paru,p1,p001, * sum,temp,zero double precision dpmpar,enorm data p1,p001,zero /1.0d-1,1.0d-3,0.0d0/ c c dwarf is the smallest positive magnitude. c dwarf = dpmpar(2) c c compute and store in x the gauss-newton direction. if the c jacobian is rank-deficient, obtain a least squares solution. c nsing = n do 10 j = 1, n wa1(j) = qtb(j) if (r(j,j) .eq. zero .and. nsing .eq. n) nsing = j - 1 if (nsing .lt. n) wa1(j) = zero 10 continue if (nsing .lt. 1) go to 50 do 40 k = 1, nsing j = nsing - k + 1 wa1(j) = wa1(j)/r(j,j) temp = wa1(j) jm1 = j - 1 if (jm1 .lt. 1) go to 30 do 20 i = 1, jm1 wa1(i) = wa1(i) - r(i,j)*temp 20 continue 30 continue 40 continue 50 continue do 60 j = 1, n l = ipvt(j) x(l) = wa1(j) 60 continue c c initialize the iteration counter. c evaluate the function at the origin, and test c for acceptance of the gauss-newton direction. c iter = 0 do 70 j = 1, n wa2(j) = diag(j)*x(j) 70 continue dxnorm = enorm(n,wa2) fp = dxnorm - delta if (fp .le. p1*delta) go to 220 c c if the jacobian is not rank deficient, the newton c step provides a lower bound, parl, for the zero of c the function. otherwise set this bound to zero. c parl = zero if (nsing .lt. n) go to 120 do 80 j = 1, n l = ipvt(j) wa1(j) = diag(l)*(wa2(l)/dxnorm) 80 continue do 110 j = 1, n sum = zero jm1 = j - 1 if (jm1 .lt. 1) go to 100 do 90 i = 1, jm1 sum = sum + r(i,j)*wa1(i) 90 continue 100 continue wa1(j) = (wa1(j) - sum)/r(j,j) 110 continue temp = enorm(n,wa1) parl = ((fp/delta)/temp)/temp 120 continue c c calculate an upper bound, paru, for the zero of the function. c do 140 j = 1, n sum = zero do 130 i = 1, j sum = sum + r(i,j)*qtb(i) 130 continue l = ipvt(j) wa1(j) = sum/diag(l) 140 continue gnorm = enorm(n,wa1) paru = gnorm/delta if (paru .eq. zero) paru = dwarf/dmin1(delta,p1) c c if the input par lies outside of the interval (parl,paru), c set par to the closer endpoint. c par = dmax1(par,parl) par = dmin1(par,paru) if (par .eq. zero) par = gnorm/dxnorm c c beginning of an iteration. c 150 continue iter = iter + 1 c c evaluate the function at the current value of par. c if (par .eq. zero) par = dmax1(dwarf,p001*paru) temp = dsqrt(par) do 160 j = 1, n wa1(j) = temp*diag(j) 160 continue call qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2) do 170 j = 1, n wa2(j) = diag(j)*x(j) 170 continue dxnorm = enorm(n,wa2) temp = fp fp = dxnorm - delta c c if the function is small enough, accept the current value c of par. also test for the exceptional cases where parl c is zero or the number of iterations has reached 10. c if (dabs(fp) .le. p1*delta * .or. parl .eq. zero .and. fp .le. temp * .and. temp .lt. zero .or. iter .eq. 10) go to 220 c c compute the newton correction. c do 180 j = 1, n l = ipvt(j) wa1(j) = diag(l)*(wa2(l)/dxnorm) 180 continue do 210 j = 1, n wa1(j) = wa1(j)/sdiag(j) temp = wa1(j) jp1 = j + 1 if (n .lt. jp1) go to 200 do 190 i = jp1, n wa1(i) = wa1(i) - r(i,j)*temp 190 continue 200 continue 210 continue temp = enorm(n,wa1) parc = ((fp/delta)/temp)/temp c c depending on the sign of the function, update parl or paru. c if (fp .gt. zero) parl = dmax1(parl,par) if (fp .lt. zero) paru = dmin1(paru,par) c c compute an improved estimate for par. c par = dmax1(parl,par+parc) c c end of an iteration. c go to 150 220 continue c c termination. c if (iter .eq. 0) par = zero return c c end C*********************************************************** subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa) integer m,n,lda,lipvt integer ipvt(lipvt) logical pivot double precision a(lda,n),rdiag(n),acnorm(n),wa(n) c ********** integer i,j,jp1,k,kmax,minmn double precision ajnorm,epsmch,one,p05,sum,temp,zero double precision dpmpar,enorm data one,p05,zero /1.0d0,5.0d-2,0.0d0/ c c epsmch is the machine precision. c epsmch = dpmpar(1) c c compute the initial column norms and initialize several arrays. c do 10 j = 1, n acnorm(j) = enorm(m,a(1,j)) rdiag(j) = acnorm(j) wa(j) = rdiag(j) if (pivot) ipvt(j) = j 10 continue c c reduce a to r with householder transformations. c minmn = min0(m,n) do 110 j = 1, minmn if (.not.pivot) go to 40 c c bring the column of largest norm into the pivot position. c kmax = j do 20 k = j, n if (rdiag(k) .gt. rdiag(kmax)) kmax = k 20 continue if (kmax .eq. j) go to 40 do 30 i = 1, m temp = a(i,j) a(i,j) = a(i,kmax) a(i,kmax) = temp 30 continue rdiag(kmax) = rdiag(j) wa(kmax) = wa(j) k = ipvt(j) ipvt(j) = ipvt(kmax) ipvt(kmax) = k 40 continue c c compute the householder transformation to reduce the c j-th column of a to a multiple of the j-th unit vector. c ajnorm = enorm(m-j+1,a(j,j)) if (ajnorm .eq. zero) go to 100 if (a(j,j) .lt. zero) ajnorm = -ajnorm do 50 i = j, m a(i,j) = a(i,j)/ajnorm 50 continue a(j,j) = a(j,j) + one c c apply the transformation to the remaining columns c and update the norms. c jp1 = j + 1 if (n .lt. jp1) go to 100 do 90 k = jp1, n sum = zero do 60 i = j, m sum = sum + a(i,j)*a(i,k) 60 continue temp = sum/a(j,j) do 70 i = j, m a(i,k) = a(i,k) - temp*a(i,j) 70 continue if (.not.pivot .or. rdiag(k) .eq. zero) go to 80 temp = a(j,k)/rdiag(k) rdiag(k) = rdiag(k)*dsqrt(dmax1(zero,one-temp**2)) if (p05*(rdiag(k)/wa(k))**2 .gt. epsmch) go to 80 rdiag(k) = enorm(m-j,a(jp1,k)) wa(k) = rdiag(k) 80 continue 90 continue 100 continue rdiag(j) = -ajnorm 110 continue return c c end C*********************************************************** subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa) integer n,ldr integer ipvt(n) double precision r(ldr,n),diag(n),qtb(n),x(n),sdiag(n),wa(n) c ********** integer i,j,jp1,k,kp1,l,nsing double precision cos,cotan,p5,p25,qtbpj,sin,sum,tan,temp,zero data p5,p25,zero /5.0d-1,2.5d-1,0.0d0/ c c copy r and (q transpose)*b to preserve input and initialize s. c in particular, save the diagonal elements of r in x. c do 20 j = 1, n do 10 i = j, n r(i,j) = r(j,i) 10 continue x(j) = r(j,j) wa(j) = qtb(j) 20 continue c c eliminate the diagonal matrix d using a givens rotation. c do 100 j = 1, n c c prepare the row of d to be eliminated, locating the c diagonal element using p from the qr factorization. c l = ipvt(j) if (diag(l) .eq. zero) go to 90 do 30 k = j, n sdiag(k) = zero 30 continue sdiag(j) = diag(l) c c the transformations to eliminate the row of d c modify only a single element of (q transpose)*b c beyond the first n, which is initially zero. c qtbpj = zero do 80 k = j, n c c determine a givens rotation which eliminates the c appropriate element in the current row of d. c if (sdiag(k) .eq. zero) go to 70 if (dabs(r(k,k)) .ge. dabs(sdiag(k))) go to 40 cotan = r(k,k)/sdiag(k) sin = p5/dsqrt(p25+p25*cotan**2) cos = sin*cotan go to 50 40 continue tan = sdiag(k)/r(k,k) cos = p5/dsqrt(p25+p25*tan**2) sin = cos*tan 50 continue c c compute the modified diagonal element of r and c the modified element of ((q transpose)*b,0). c r(k,k) = cos*r(k,k) + sin*sdiag(k) temp = cos*wa(k) + sin*qtbpj qtbpj = -sin*wa(k) + cos*qtbpj wa(k) = temp c c accumulate the tranformation in the row of s. c kp1 = k + 1 if (n .lt. kp1) go to 70 do 60 i = kp1, n temp = cos*r(i,k) + sin*sdiag(i) sdiag(i) = -sin*r(i,k) + cos*sdiag(i) r(i,k) = temp 60 continue 70 continue 80 continue 90 continue c c store the diagonal element of s and restore c the corresponding diagonal element of r. c sdiag(j) = r(j,j) r(j,j) = x(j) 100 continue c c solve the triangular system for z. if the system is c singular, then obtain a least squares solution. c nsing = n do 110 j = 1, n if (sdiag(j) .eq. zero .and. nsing .eq. n) nsing = j - 1 if (nsing .lt. n) wa(j) = zero 110 continue if (nsing .lt. 1) go to 150 do 140 k = 1, nsing j = nsing - k + 1 sum = zero jp1 = j + 1 if (nsing .lt. jp1) go to 130 do 120 i = jp1, nsing sum = sum + r(i,j)*wa(i) 120 continue 130 continue wa(j) = (wa(j) - sum)/sdiag(j) 140 continue 150 continue c c permute the components of z back to components of x. c do 160 j = 1, n l = ipvt(j) x(l) = wa(j) 160 continue return c c end C***********************************************************