FLG - a fitfun program example accessing COMMON variables

R. Ghosh, October 1999

FLG fits different model peak functions to a spectrum. It shows how parameters can be used to control choice of model, as an alternative to being a fitting variable. Two overlapping peaks are treated here, but it is simple to extend the model to include additional peaks, and add convolutions with resolution functions etc.

This example shows how a production program can be developed from a prototype, using the additional J command. The programmer supplies the routines


FTPXTP        to do more elegant plotting of final results etc.
FTXHLP        to give addditional help information

MAIN Program and extra help information
      PROGRAM FLG
C****** fitting diffraction data...
C f77 -o flg flg.f difin.f cald.f caldx.f libfitfun.a libpgplot.a
C -lX11
      EXTERNAL DIFIN,CALD
C     DIFIN reads in a diffraction pattern
C     CALD  calculates sum of up to two peaks
      COMMON/WORK/W(160000)
      COMMON/TITLES/NAMES(40),TX,TY
      COMMON/TITLEP/NPARAS
      COMMON/VERSION/VERP
      CHARACTER*8 NAMES
      CHARACTER*4 PNAM
      CHARACTER*20 TX,TY
      DATA NAMES/ 'FLAT BGD','slope BG',
     1'P1 TYPE ','P1 AREA ','P1 WIDTH','P1 POSN ',
     1'P2 TYPE ','P2 AREA ','P2 WIDTH','P2 POSN ',
     130*'        '/
      DATA TX /'Angle'/
      DATA TY /'Counts'/
      NPARAS=10
      VERP=1.1
      WRITE(6,1) VERP
1     FORMAT(' FLG - fitting diffraction data...'/
     1/' VERSION ',F4.1,'  -  R.E. GHOSH  October ''99'//
     1' PROGRAM FITS SINGLE peaks to Counts'//
     1' PEAK TYPE :'/
     1'           1 SINGLE LORENTZIAN'/
     1'           2 SINGLE GAUSSIAN'/
     1//' COMMAND J FOR COMPONENT  PLOTS'/)
C***** LIMIT FOR TESTS
      PNAM='flgf'
C***** OPEN DATA FILES
      CALL FTFUNS(PNAM,DIFIN,CALD)
      END
      SUBROUTINE FTXHLP(IOTTY)
      WRITE(IOTTY,1)
1     FORMAT(/' FLG'
     1' Program fits peaks to diffraction data'//
     1' PEAK TYPE 1 SINGLE LORENTZIAN'/
     1'           2 SINGLE GAUSSIAN'/
     1/' COMMAND J FOR COMPONENT  PLOTS'/)
      RETURN
      END
Data Reading routine

      SUBROUTINE DIFIN(NPNTS,XOBS,YOBS,YERR,NTEXT)
C***** Reads in a diffraction scan
      DIMENSION XOBS(1),YOBS(1),YERR(1)
      COMMON/INDIFC/FNAME
C***** Save - may be useful for titles
      CHARACTER*50 NTEXT
      CHARACTER*20 FNAME
10    WRITE(6,1)
1     FORMAT(//'   Give input filename : ',$)
      READ(5,2,END=100) FNAME
C***** exit with control-z is possible here though plots are
C      NOT TERMINATED
2     FORMAT(A)
      OPEN(UNIT=10,FILE=FNAME,ACCESS='SEQUENTIAL',STATUS='OLD',ERR=95)
C***** note unit is not 1
C***** simply formatted file - skip headers
      READ(10,3,ERR=99,END=99) NTEXT
      READ(10,4,ERR=99,END=99)
      READ(10,4,ERR=99,END=99)
      READ(10,4,ERR=99,END=99)
4     FORMAT(1x)
3     FORMAT(A)
C***** this can serve as a title for the spectrum
      NPNTS=0
      DO 5 I=1,1000
      READ(10,6,ERR=99,END=99) XOBS(I),YOBS(I)
      yerr(i)=sqrt(yobs(i))
C***** set errors to sqrt(counts)
6     FORMAT(3G)
      NPNTS=I
5     CONTINUE
99    CLOSE(UNIT=10)
      IF(NPNTS.LE.5) GO TO 12
      WRITE(6,9) NTEXT
9     FORMAT(' TITLE: ',a/)
      WRITE(6,11) NPNTS
11    FORMAT(1X,I4,' data have been read in ')
      RETURN
12    WRITE(6,98)
98    FORMAT(' INSUFFICIENT DATA TO CONTINUE')
C***** do not return until some good data are read
      GO TO 10
95    close(10)
      write(6,94)
94    format(' File not found ')
      go to 10
100   STOP 'END OF TERMINAL INPUT'
      END
Calculation routine

Note that the component curves are conserved in COMMON storage.


      SUBROUTINE CALPD(NP,P,NFIT,ANG,YUSE,YRUSE,YCALC,F)
C***** Calculates sum of two peaks for fitfun fit to diff data
C     INPUT PARAMETERS
C     P 1       Flat background
C     P 2       slope        
C     P 3       PEAK TYPE (1=LOR,2=GAU)
C     P 4       PEAK HEIGHT
C     P 5       PEAK WIDTH
C     P 6       PEAK POSN
C     P 7       PEAK TYPE FOR 2ND COMPONENT... ETC...
C***** simplified example program derived from FIRR and TAPPIT
      COMMON/EXTPC/YY(250,3),YYIN(2),INPK
C***** Save calculated curves for additional treatment...
      DIMENSION P(NP),ANG(NFIT),YUSE(NFIT),YRUSE(NFIT),YCALC(NFIT)
     1,F(NFIT)
      DIMENSION Y(250)
C
C***** Calculate y - no resolution effects
C
      ALN2=ALOG(2.0)
      SLN2=SQRT(ALN2)
      PI=3.141592
      SPI=SQRT(PI)
      YYIN(1)=0.
      YYIN(2)=0.
      DO 1000 I=1,NFIT
      YY(I,1)=0.
      YY(I,2)=0.
      YY(I,3)=0.
      YCALC(I)=0.
1000  CONTINUE
1     CONTINUE
      INPK=0
      DO 2 III=3,7,4
C***** FOR ANY PEAKS PRESENT
      IF(P(III).EQ.0) GO TO 2
C***** NOT THIS ONE
      INPK=INPK+1
      YYIN(INPK)=1.
      IH=III+1
C***** Area
      IW=III+2
C***** Width
      IP=III+3
C***** Position
      J=P(III)
      IF(J.LT.1.OR.J.GT.2) GO TO 2
C***** Other types of peak not treated
      GO TO (20,30) J
20    CONTINUE
C***** Lorentzian function
      DO 21 I=1,NFIT
      Y(I)=0.
      ANGF=ANG(I)
      FF=(ANGF-P(IP))**2
      Y(I)=Y(I)+P(IH)*P(IW)/(FF+P(IW)**2)/PI
21    CONTINUE
      GO TO 202
30    CONTINUE
C***** Gaussian function
      DO 31 I=1,NFIT
      Y(I)=0.
      ANGF=ANG(I)
      FF=(ANGF-P(IP))**2
      Y(I)=Y(I)+P(IH)*EXP(-FF*ALN2/P(IW)**2)*SLN2/(SPI*P(IW))
31    CONTINUE
      GO TO 202
202   DO 204 I=1,NFIT
      YY(I,INPK)=Y(I)
C***** Each component is saved for possible later use
204   CONTINUE
2     CONTINUE
      DO 170 I=1,NFIT
      YY(I,3)=P(1)+P(2)*ANG(I)
C***** Flat background+slope
      YCALC(I)=YY(I,1)+YY(I,2)+YY(I,3)
170   CONTINUE
      DO 200 I=1,NFIT
      FFFF=YUSE(I)-YCALC(I)
      F(I)=FFFF/YRUSE(I)
200   CONTINUE
      RETURN
      END

Additional output can be provided using the J command and the routine FTEXTP:
      SUBROUTINE FTEXTP
      PARAMETER (MAXDAT=21000)
      PARAMETER (MAXPAR=40)
C***** EXTERNAL SUBROUTINE TO PLOT Component RESULTS
      COMMON/PLTLIM/XMIN,XMAX,YMIN,YMAX,XSCALE(2),YSCALE(2)
      COMMON/PLTLIC/NTX,NTY,NTEX
      COMMON/FUNCTN/NPAR,PARM(MAXPAR),STEP(MAXPAR),SCALE(MAXPAR)
      COMMON/OBSDAT/NPNT,XUSE(MAXDAT),YUSE(MAXDAT),YUSR(MAXDAT)
      COMMON/CALDAT/YCALC(MAXDAT),F(MAXDAT),IITER
      COMMON/REDITC/XOBS(MAXDAT),YOBS(MAXDAT),YROBS(MAXDAT),YER(MAXDAT)
     1,NDIN
      COMMON/INDIFC/FNAME
      COMMON/TITLES/NAMES(MAXPAR),TX,TY
      COMMON/CHICHI/RF,RFOLD,CHI,CHIOLD,PROB,NFREE,NFITNO
      COMMON/TITLEX/EXTRAT
      COMMON/EXTPC/YY(250,3),YYIN(2),INPK
c***** other commons with useful data are
      COMMON/EXTPLT/XPL(200),YPL(200),WERR(200),WRK(200),WK1(200),NEXTRA
      COMMON/FITLIM/XLIMS(6),NSET,MINF,MAXF,NFIT
c***** which have the x range limits "O", and extrapolated values
      CHARACTER*8 NAMES
      CHARACTER*20 TX,TY
      CHARACTER*20 FNAME
      CHARACTER NTX*20,NTY*20,NTEX*50,NTDUM*50,EXTRAT*50
      CHARACTER*1 ANS
      CHARACTER PNAM*4
      CHARACTER*80 PLTXT(2),BTXT
      CHARACTER*4 PTYP(2),ACTTYP
      CHARACTER*6 COLS(2)
      INTEGER MAPCOL(2)
      REAL YLOC(MAXDAT),FLOC(MAXDAT),YDIF(250)
      DATA MAPCOL/2,3/
      DATA PTYP/'LOR_','GAUS'/
      DATA COLS/'RED   ','GREEN '/
      PLTXT(1)=' '
      PLTXT(2)=' '
      BTXT=' '
      NTDUM=' '
C***** perform calculation for all points once again
      CALL CALD(NPAR,PARM,NDIN,XOBS,YOBS,YROBS,YCALC,F)
      OFFSET=0.8*(YSCALE(2)-YSCALE(1))+YSCALE(1)
C***** Difference 
      DO 6000 I=1,NPNT
      YDIF(I)=OFFSET+YCALC(I)-YOBS(I)
      YY(I,1)=YY(I,1)+YY(I,3)
      YY(I,2)=YY(I,2)+YY(I,3)
C***** PLot components on top of notional background
6000  CONTINUE
C***** Load up titles
      LIN=0
      IPARAT=-1
      DO 7000 I=1,2
      IPARAT=IPARAT+4
      IF(YYIN(I).EQ.0.) GO TO 7000
C***** Component absent
      LIN=LIN+1
      ITYPX=PARM(IPARAT)
      ACTTYP=' '
      IF(ITYPX.GT.0.AND.ITYPX.LE.2) ACTTYP=PTYP(ITYPX)
      WRITE(PLTXT(LIN),7001) ACTTYP,
     1NAMES(IPARAT+1),PARM(IPARAT+1),
     1NAMES(IPARAT+2),PARM(IPARAT+2),
     1NAMES(IPARAT+3),PARM(IPARAT+3),COLS(I)
7001  FORMAT(A4,3(1X,A8,F10.5),2X,A6)
7000  CONTINUE
      IF(CHI.NE.0.) WRITE(BTXT,7010)NFITNO,CHI
7010  FORMAT(' FIT #',I3,40X,'CHISQUARED ',F8.2)
C***** Graph out
C***** Frame and text
      IER=0
      CALL RSPLT(XSCALE,YSCALE,YSCALE,2,0,IER,NTX,NTY,NTDUM)
C***** Data points
      CALL RSPLT(XOBS,YOBS,YROBS,NDIN,-2,IER,NTX,NTY,NTEX)
C***** Fit
      CALL RSPLT(XOBS,YCALC,YCALC,NDIN,-100,IER,NTX,NTY,NTEX)
C***** Difference
      CALL RSPLT(XOBS,YDIF,YDIF,NDIN,-100,IER,NTX,NTY,NTEX)
C***** First component red
      CALL PGQCI(IOLD)
C***** Save initial colours
      DO 7500 I=1,2
      IF(YYIN(I).EQ.0.) GO TO 7500
      CALL PGSCI(MAPCOL(I))
      IF(YYIN(1).NE.0)
     1CALL RSPLT(XOBS,YY(1,I),YY(1,I),NDIN,-100,
     1IER,NTX,NTY,NTEX)
7500  CONTINUE
      CALL PGSCI(IOLD)
C***** Background
      CALL RSPLT(XOBS,YY(1,3),YY(1,3)
     1,NDIN,-200,IER,NTX,NTY,NTEX)
c***** titles
      CALL PGMTXT('T',7.,0.,0.,FNAME//NTEX)
      CALL PGMTXT('T',5.,0.,0.,PLTXT(1))      
      CALL PGMTXT('T',4.,0.,0.,PLTXT(2))
      CALL PGMTXT('T',3.,0.,0.,EXTRAT)      
      CALL PGMTXT('B',2.5,0.5,0.5,NTX)
      CALL PGMTXT('T',.75,0.5,0.5,BTXT)
      WRITE(6,300)
      CALL PGSCI(IOLD)
300   FORMAT(' PostScript Output (Y/N) : ',$)
      READ(5,301,END=302) ANS
301   FORMAT(A1)
302   CONTINUE
      IF(ANS.EQ.'Y'.OR.ANS.EQ.'y') GO TO 303
      RETURN
303   CONTINUE
      CALL PGSLCT(2,LAST,IER)
C***** Graph out
C***** FRAME AND TEXT
      IER=0
      CALL RSPLT(XSCALE,YSCALE,YSCALE,2,0,IER,NTX,NTY,NTDUM)
C***** Data points
      CALL RSPLT(XOBS,YOBS,YROBS,NDIN,-2,IER,NTX,NTY,NTEX)
C***** Fit
      CALL RSPLT(XOBS,YCALC,YCALC,NDIN,-100,IER,NTX,NTY,NTEX)
C***** Difference
      CALL RSPLT(XOBS,YDIF,YDIF,NDIN,-100,IER,NTX,NTY,NTEX)
      CALL PGQCI(IOLD)
      DO 7600 I=1,2
      IF(YYIN(I).EQ.0.) GO TO 7600
      CALL PGSCI(MAPCOL(I))
      IF(YYIN(1).NE.0)
     1CALL RSPLT(XOBS,YY(1,I),YY(1,I),NFIT,-100,
     1IER,NTX,NTY,NTEX)
7600  CONTINUE
      CALL PGSCI(IOLD)
C***** Background
      CALL RSPLT(XOBS,YY(1,3),YY(1,3)
     1,NDIN,-200,IER,NTX,NTY,NTEX)
C***** Titles
      CALL PGMTXT('T',7.,0.,0.,FNAME//NTEX)
      CALL PGMTXT('T',5.,0.,0.,PLTXT(1))      
      CALL PGMTXT('T',4.,0.,0.,PLTXT(2))      
      CALL PGMTXT('T',3.,0.,0.,EXTRAT)      
      CALL PGMTXT('B',2.5,0.5,0.5,NTX)
      CALL PGMTXT('T',.75,0.5,0.5,BTXT)
      CALL PGSLCT(LAST,LLAST,IER)
      RETURN
      END
The graphics library routines beginning with the letters PG are described in the PGPLOT library page . RSPLT in the fitfun library is identical to that in the librlib.a library.

A test dataset, follows together with an example of the program...

file: test
  1.300  1.500  0.500 to   1.500  1.300  0.500
 02321.dif
 20000.0 0.167
 2
    1.300000       153.0000    
    1.302500       144.0000    
    1.305000       176.0000    
    1.307500       175.0000    
    1.310000       173.0000    
    1.312500       166.0000    
    1.315000       179.0000    
    1.317500       201.0000    
    1.320000       195.0000    
    1.322500       182.0000    
    1.325000       178.0000    
    1.327500       162.0000    
    1.330000       173.0000    
    1.332500       189.0000    
    1.335000       188.0000    
    1.337500       171.0000    
    1.340000       168.0000    
    1.342500       210.0000    
    1.345000       183.0000    
    1.347500       177.0000    
    1.350000       218.0000    
    1.352500       194.0000    
    1.355000       210.0000    
    1.357500       208.0000    
    1.360000       219.0000    
    1.362500       198.0000    
    1.365000       220.0000    
    1.367500       243.0000    
    1.370000       241.0000    
    1.372500       250.0000    
    1.375000       270.0000    
    1.377500       281.0000    
    1.380000       293.0000    
    1.382500       281.0000    
    1.385000       310.0000    
    1.387500       323.0000    
    1.390000       331.0000    
    1.392500       370.0000    
    1.395000       454.0000    
    1.397500       548.0000    
    1.400000       704.0000    
    1.402500       626.0000    
    1.405000       494.0000    
    1.407500       416.0000    
    1.410000       424.0000    
    1.412500       375.0000    
    1.415000       376.0000    
    1.417500       396.0000    
    1.420000       356.0000    
    1.422500       385.0000    
    1.425000       365.0000    
    1.427500       317.0000    
    1.430000       326.0000    
    1.432500       354.0000    
    1.435000       286.0000    
    1.437500       300.0000    
    1.440000       298.0000    
    1.442500       285.0000    
    1.445000       266.0000    
    1.447500       270.0000    
    1.450000       288.0000    
    1.452500       248.0000    
    1.455000       263.0000    
    1.457500       210.0000    
    1.460000       235.0000    
    1.462500       231.0000    
    1.465000       230.0000    
    1.467500       238.0000    
    1.470000       208.0000    
    1.472500       234.0000    
    1.475000       219.0000    
    1.477500       226.0000    
    1.480000       232.0000    
    1.482500       233.0000    
    1.485000       220.0000    
    1.487500       192.0000    
    1.490000       212.0000    
    1.492500       233.0000    
    1.495000       214.0000    
    1.497500       205.0000    
    1.500000       218.0000    

The script for the test session follows:

% flg

 FLG - fitting diffraction data...

 VERSION  1.1  -  R.E. GHOSH  October '99

 PROGRAM FITS SINGLE peaks to Counts

 PEAK TYPE :
           1 SINGLE LORENTZIAN
           2 SINGLE GAUSSIAN


 COMMAND J FOR COMPONENT  PLOTS

 flgf               version  1.1
 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
 flgf>r


   Give input filename : test
 TITLE:   1.300  1.500  0.500 to   1.500  1.300  0.500    

   81 data have been read in 
 THERE ARE  81 DATA POINTS IN CURRENT FITTING RANGE
 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
 flgf>v
 FITFUN 5.4  TITLE:  1.300  1.500  0.500 to   1.500  1.300  0.500    
 FITTING Y : Counts                 VERSUS X : Angle               
 NUMBER  PARAMETER   VALUE    (  OLD VALUE )     STEP     % DEVIATION 
    1   FLAT BGD     180.0    (   180.0    )   1.000           0.0
    2   slope BG     10.00    (   10.00    )   1.000           0.0
    3   P1 TYPE      2.000    (   2.000    )  0.0000E+00       0.0
    4   P1 AREA      2.500    (   2.500    )   1.000           0.0
    5   P1 WIDTH    0.4000E-02(  0.4000E-02)   1.000           0.0
    6   P1 POSN      1.400    (   1.400    )  0.5000E-01       0.0
    7   P2 TYPE      1.000    (   1.000    )  0.0000E+00       0.0
    8   P2 AREA      20.00    (   20.00    )   1.000           0.0
    9   P2 WIDTH    0.3000E-01(  0.3000E-01)   1.000           0.0
   10   P2 POSN      1.400    (   1.400    )  0.5000E-01       0.0
  100     MAXIMUM STEP    100.0       300  SUBSTEP   0.10000
  200     ACCURACY      0.1000E-01
 THERE ARE  81 POINTS IN THE CURRENT RANGE SET BY "ONLY" COMMAND.
 CURRENT LIMITS (IF RANGE IS NON-ZERO) ARE : 
  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00
 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
 flgf>f

 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
flgf>f 100
 FITTING.....

 .....ENDED
 

 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
 flgf>j

 PostScript Output (Y/N) : y
 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
 flgf>l
 FITFUN 5.4  TITLE:   1.300  1.500  0.500 to   1.500  1.300  0.500    
                                                                      
 FITTING Y : Counts                         FIT NUMBER    1
   VERSUS X : Angle               
 NUMBER  PARAMETER   VALUE    (  OLD VALUE )    STEP      % DEVIATION 
    1   FLAT BGD    -148.3    (   180.0    )   1.000         -24.6
    2   slope BG     223.7    (   10.00    )   1.000         325.1
    3   P1 TYPE      2.000    (   2.000    )  0.0000E+00       0.0
    4   P1 AREA      2.555    (   2.500    )   1.000           8.8
    5   P1 WIDTH    0.3934E-02(  0.4000E-02)   1.000           9.5
    6   P1 POSN      1.400    (   1.400    )  0.5000E-01       0.0
    7   P2 TYPE      1.000    (   1.000    )  0.0000E+00       0.0
    8   P2 AREA      22.11    (   20.00    )   1.000           4.2
    9   P2 WIDTH    0.3075E-01(  0.3000E-01)   1.000           3.3
   10   P2 POSN      1.409    (   1.400    )  0.5000E-01       0.1
  100     MAXIMUM STEP    100.0     SUBSTEP    0.10000
  200     ACCURACY       0.1000E-01
 CURRENT RESIDUAL VALUE :  0.10291    
 THERE ARE  81 POINTS IN THE CURRENT RANGE SET BY "ONLY" COMMAND.
 CURRENT LIMITS (IF RANGE IS NON-ZERO) ARE : 
  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00  0.000E+00 TO 0.000E+00


 Mean R-factor is now   4.9    %     (old value   0.00E+00%)
 Reduced CHI-squared is now   0.95     for    73 degrees of freedom
 (old value   0.00E+00)
 Chi-squared probability is    0.595    

 TYPE HELP OR OPTION: H,R,D,X,Y,V,F,O,P,L,W,C,E,J,S,T,Z
 flgf>e
 DO YOU WISH TO SAVE THE CURRENT PARAMETER VALUES ? (Y)
n
 Closing listing file flgf029.lis
 %PGPLOT, Completing PostScript file flgf028.ps
is1 16% lp flgf029.lis
request id is dj_rg-792 (1 file(s))
is1 17% lp flgf028.ps
request id is dj_rg-793 (1 file(s))