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 informationMAIN Program and extra help informationPROGRAM 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' call prelude(pnam) c***** prepares for use with Clickfit-GUI 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 ENDData Reading routineSUBROUTINE 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' ENDCalculation routineNote 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 ENDAdditional 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 ENDThe 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.0000The 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>jPostScript 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))