      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,'      <I(R40)>:',oI40
      WRITE(*,'(A,F10.3,A,F10.3)')'  R50:',oR50,'      <I(R50)>:',oI50
      WRITE(*,'(A,F10.3,A,F10.3)')'  R60:',oR60,'      <I(R60)>:',oI60
      WRITE(*,'(A,F10.3,A,F10.3)')'  R70:',oR70,'      <I(R80)>:',oI70
      WRITE(*,'(A,F10.3,A,F10.3)')'  R80:',oR80,'      <I(R80)>:',oI80
      WRITE(*,'(A,F10.3,A,F10.3)')'  R90:',oR90,'      <I(R90)>:',oI90
      WRITE(*,'(A,F10.3,A,F10.3)')' R100:',oR100,'     <I(R100)>:',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),'      <I(R40)>:',Ix(4)
      WRITE(*,'(A,F10.3,A,F10.3)')'  R50:',Rx(5),'      <I(R50)>:',Ix(5)
      WRITE(*,'(A,F10.3,A,F10.3)')'  R60:',Rx(6),'      <I(R60)>:',Ix(6)
      WRITE(*,'(A,F10.3,A,F10.3)')'  R70:',Rx(7),'      <I(R80)>:',Ix(7)
      WRITE(*,'(A,F10.3,A,F10.3)')'  R80:',Rx(8),'      <I(R80)>:',Ix(8)
      WRITE(*,'(A,F10.3,A,F10.3)')'  R90:',Rx(9),'      <I(R90)>:',Ix(9)
      WRITE(*,'(A,F10.3,A,F10.3)')' R100:',Rx(10),
     &'     <I(R100)>:',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&amp;data_type=PDF_HIGH&amp;whole_paper=YES&amp;type=PRINTER&amp;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***********************************************************      
