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=2100)
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
Legal notice