iraf.net

One of the long-forgotten ADASS newsgroups was adass.iraf.sources, the point of which was to allow users to post small bits of code others may find useful. While we work to decide how we'll organize something similar we'll simply post the code in a new topic.

Weidong Li at UCBerkeley asked recently how one can use the coefficients produced by CURFIT to evaluate the fit at arbitrary points. After a quick pointer, and in what we hope will be an iraf.net tradition, he coded a solution and agreed to share it here. Thanks Weidong!

CC      Example program to evaluate curfit Spline3 function
C       at any arbitrary x value (x must be between xmax and xmin)
C       
C       To use the program:
C
C       (1) the spline3 order should be in a file
C           named coeff.dat. 
C       (3) The X values at which the Spline3 function should be 
C           evaluated should be in a file named input.dat. Make sure
C           all the X values are between xmax and xmin! The maximum
C           number of X values is 10000
C       (4) The program will ask for the following parameters:
C           the Spline3 order; the minimun X value used in the 
C           fit; the maximum X value used in the fit
C       (5) the output is in a file called output.dat
C       
C       Written by Weidong Li at UC Berkeley 2005-Dec-05
C       weidong@astron.berkeley.edu
C       



        program evaluate_a_spline3

        real coeff(100),basis(4)

        open(21,file='coeff.dat')
        open(22,file='input.dat')
        open(23,file='output.dat')

        write(*,*)
        write(*,*) 'Input the order of the fit:'
        read(*,*) np
        write(*,*) 'Input the mimimum X used in the fit: '
        read(*,*) xmin
        write(*,*) 'Input the maximum X used in the fit: '
        read(*,*) xmax

C       There must be np+3 coefficients in the file coeff.dat
C       otherwise the program will quit!

        ncoeff=np+3

        do i=1,ncoeff
           read(21,*) coeff(i)
        enddo

        dx  = np/(xmax-xmin)

       do n=1,10000
        read(22,*,end=999) x
        sx=(x-xmin)*dx
        ileft=min(int(sx),np)

        sx = max (0, min(1.0,sx-ileft))
        tx = max (0, min(1.0,1.0-sx))

        basis(1) = tx*tx*tx
        basis(2) = 1 + tx * ( 3 + tx*(3-3*tx))
        basis(3) = 1 + sx * ( 3 + sx*(3-3*sx))
        basis(4) = sx*sx*sx

        yfit=0.0
        do j=ileft+1,ileft+4
           yfit=yfit+coeff(j)*basis(j-ileft)
        enddo

        write(23,1001) x,yfit
       enddo

999    continue

1001   format(5(3x,f9.4))

        stop
        end