Close This Window

Please download official ILL logos here

 

For using on the web or on a screenFor printing in high resolutionWhite version, for dark backgrounds

Download PNG

Download AI

Download white PNG

Download JPG

 

Download white AI

FLG - a fitfun program example

The Computing for Science (CS) group supports ILL scientists, students and visitors in a number of activities including data analysis, instrument simulation and sample simulation.

Back to ILL Homepage
English French Deutsch 

All Software

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=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
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))


To top