to run this program, compile progrw.f and link to make progrw.exeg77 -o progrw progrw.f and run progrw 10 lines have been read line 1 has 1.6860990E-03 464.5000 10.77613 line 2 has 4.1995151E-03 493.0833 6.410170 line 3 has 6.5153618E-03 476.0625 5.454714 line 4 has 9.1368612E-03 444.2500 3.983224 line 5 has 1.1507930E-02 496.5000 4.210955 line 6 has 1.3852740E-02 512.4167 3.772770 line 7 has 1.6353959E-02 510.4375 3.260999 line 8 has 1.8697230E-02 503.2954 3.382090 line 9 has 2.1087030E-02 493.0333 2.866570 line 10 has 2.3404090E-02 507.0179 3.008968Clearly one wants more flexibility than just reading the file "mydata.dat". In the next example the reading and writing activities are localised in subroutines. Because there might be errors in the data it is useful to have a way of showing this.
Second example using reading and writing subroutinesProgram progrwm.f - written with Notepad program progrwm c***** a simple program to read and write data using subroutines c R. Ghosh, ILL, version 1 October 2004 c g77 -o progrwm progrwm.f daread.f dawrit.f parameter (inmax=100) c***** maximum number of lines to be read from file real x(inmax),y(inmax),yer(inmax) c***** define arrays of variables x,y,yer call daread(x,y,yer,in,inmax) c***** routine reads in from file iout=6 call dawrit(x,y,yer,in,iout) c***** writes out to unit 6 (terminal) end and a read routine daread.f subroutine daread(xx,yy,yr,in,inmax) c****** routine to read from a 3 column data file c xx out first column value c yy out second column value c yr out third column value c in out number of lines read c inmax in maximum lines expected c R. Ghosh, Ill, October 2004 real xx(inmax),yy(inmax),yr(inmax) common/infile/fname character*100 fname c***** the filename to be read has a maximum length of 100 characters write(6,1) 1 format(' Give input filename : ',$) read(5,2) fname 2 format(a) open(unit=2,file=fname) do i=1,inmax read(2,*,err=25,end=20) xx(i),yy(i),yr(i) j=i end do 20 continue close(unit=2) c***** not required further - close input write(6,3) j 3 format(i3,' lines have been read in') in=j return 25 write(6,26)j 26 format(' daread - error after line ',i4) c***** when reading data it is useful to offer guidance on c potential errors therein stop end and a write routine dawrit.f subroutine dawrit(xx,yy,yr,nin,nunit) c****** routine to write out a 3 column data file c if nunit is 6 then output is to terminal c otherwise output is to a named file c c xx in first column value c yy in second column value c yr in third column value c nin in number of lines <=100 c nunit in output unit c R. Ghosh, ILL, October 2004 real xx(nin),yy(nin),yr(nin) common/oufile/fname character*120 fname c***** the filename to be read has a maximum length of 120 characters if(nin.gt.100) stop ' dawrit- cannot write more than 100 data' c***** all files treated by these routines should have less than 101 lines if(nunit.ne.6) then write(6,1) 1 format(' Give output filename : ',$) read(5,2) fname 2 format(a) end if c***** open a new file if(nunit.ne.6) open(unit=nunit,file=fname,status='new') c***** could use 'unknown' rather than 'new' if allow overwriting do i=1,nin write(nunit,*) xx(i),yy(i),yr(i) j=i end do if(nunit.ne.6) close(unit=nunit) c***** finished writing file output - close file write(6,3) j 3 format(i3,' lines have been written') return end putting it all together g77 -o progrwm progrwm.f daread.f dawrit.f progrwm Give input filename : mydata.dat 10 lines have been read in 1.6860990E-03 464.5000 10.77613 4.1995151E-03 493.0833 6.410170 6.5153618E-03 476.0625 5.454714 9.1368612E-03 444.2500 3.983224 1.1507930E-02 496.5000 4.210955 1.3852740E-02 512.4167 3.772770 1.6353959E-02 510.4375 3.260999 1.8697230E-02 503.2954 3.382090 2.1087030E-02 493.0333 2.866570 2.3404090E-02 507.0179 3.008968 10 lines have been written
Subtracting two sets of data is the next exampleprogram progsub c***** a simple program to read two data files and subtract c the second from the first, using subroutines c R. Ghosh, ILL, version 1 October 2004 c g77 -o progsub progsub.f daread.o dawrit.o parameter (inmax=100) c***** maximum number of lines to be read from file real x1(inmax),y1(inmax),yer1(inmax) real x2(inmax),y2(inmax),yer2(inmax) real x3(inmax),y3(inmax),yer3(inmax) c***** define arrays of variables x,y,yer call daread(x1,y1,yer1,in1,inmax) c***** routine reads in from file call daread(x2,y2,yer2,in2,inmax) c***** trivial check for compatibility if(in2.ne.in1) stop ' input files have differing lengths!' write(6,5) 5 format(' Give multiplier for second set of data [1.0] :',$) read(5,6) fmul 6 format(g12.6) if(fmul.eq.0) fmul=1.0 c c***** calculate difference and new error values c do i=1,in1 x3(i)=x1(i) c***** take first data as reference y3(i)=y1(i)-fmul*y2(i) yysq=yer1(i)**2+(fmul*yer2(i))**2 yer3(i)=sqrt(yysq) end do c c***** output results c iout=6 call dawrit(x3,y3,yer3,in1,iout) c***** writes out to unit 6 (terminal) iout=2 c***** output into a new file call dawrit(x3,y3,yer3,in1,iout) end Reusing the routine from above: g77 -o progsub progsub.f daread.f dawrit.f or, because the compiled subroutines remain there: g77 -o prgsub progsub.f daread.o dawrit.o and test it.. progsub Give input filename : mydata.dat 10 lines have been read in Give input filename : mydata.dat 10 lines have been read in Give multiplier for second set of data [1.0] :-1.0 1.6860990E-03 929.0000 15.23975 4.1995151E-03 986.1666 9.065350 6.5153618E-03 952.1250 7.714130 9.1368612E-03 888.5000 5.633129 1.1507930E-02 993.0000 5.955190 1.3852740E-02 1024.833 5.335503 1.6353959E-02 1020.875 4.611749 1.8697230E-02 1006.591 4.782998 2.1087030E-02 986.0666 4.053942 2.3404090E-02 1014.036 4.255323 10 lines have been written Give output filename : newdata.dat 10 lines have been writtenas a simple check the y data have doubled, and the error increased by a factor of 1.414 which is correct if the data are from two uncorrelated sets of counting events(SQRT(2.0))
Plotting data using library routine RSPLT from library librlib.a which, in turn, calls Prof. Tim Pearson's PGPLOT libraryprogram daplo c example of command line driven plotting program c R. Ghosh ILL, October 2004, version 1 c SGI f77 -o daplo daplo.f daread1.f librlib.a libpgplot.a -lX11 cMinGW g77 -o daplo daplo.f daread1.f librlib.a libpggw520.a -wl,subsystem,console -mwindows c c***** parameter (inmax=100) real x(inmax),y(inmax),yer(inmax) common/infile/filnam character*100 filnam,str character*50 title noargs=iargc() c***** gives number of fields in command line if(noargs.lt.1.or.noargs.gt.6) then write(*,*) 1' plot data : usage daplo filename ymin ymax xmin xmax' call exit end if filnam=' ' str=' ' c***** set the whole string str to blanks call getargp(1,str) c***** this is first character string in command line - file name filnam=str call daread1(x,y,yer,nin,inmax) c***** daread1 will get filename from filnam iidev=1 call pgnopn(iidev,' ',ier) c***** open graphics if (ier.ne.0) stop 'no screen device PGPLOT_ILL_DEV_1' c***** put filename into title for plot title=str call rsplt(x,y,yer,nin,2,1,'X','Y',title) call pgnend(iidev) end subroutine daread1(xx,yy,yr,in,inmax) c****** routine to read from a 3 column data file c xx out first column value c yy out second column value c yr out third column value c in out number of lines read c inmax in maximum lines expected c Input filename set in common/infile/fname c R. Ghosh, Ill, October 2004 real xx(inmax),yy(inmax),yr(inmax) common/infile/fname character*100 fname c***** the filename to be read has a maximum length of 100 characters open(unit=2,file=fname) do i=1,inmax read(2,*,end=20) xx(i),yy(i),yr(i) j=i end do 20 continue close(unit=2) c***** not required further - close input write(6,3) j 3 format(i3,' lines have been read in') in=j return end build it g77 -o daplo daplo.f daread1.f librlib.a libpggw520.a -wl,subsystem,console -mwindows now we can run this directly from the command line.. daplo sasteflon.datThe batch file mingwvars.bat may be used as a shortcut to setting up a CMD program window used for development, modifying this shortcut to start in an appropriate test directory. Part of the g77 command line is setup as an environment variable %PGFLAG% simplifying the above command to:
g77 -o daplo daplo.f daread1.f librlib.a %PGFLAG%
The command for building the program is getting rather complicated and a useful tool is gmake (GNU make) which uses a text "makefile" to build only those components which have been changed. The "target" (result) is denoted by the name+colon followed by the files used to create it. The following line has the appropriate command, starting with a colon. In the MinGW distribution the make program is initially named mingw32-make.exenotepad makefile all: daplo progsub progrwm daplo.exe: daplo.f daread1.o g77 -o daplo daplo.f daread1.o librlib.a libpggw520.a -wl,subsystem,console -mwindows daread1.o: daread1.f g77 -c daread1.f progsub.exe: progsub.f daread.o dawrit.o g77 -o progsub progsub.f daread.o dawrit.o progrw.exe: progrw.f daread.o dawrit.o g77 -o progrw progrw.f daread.o dawrit.o daread.o: daread.f g77 -c daread.f dawrit.o: dawrit.f f77 -c dawrit.f clean: del *.o The command make daplo makes daplo.exe The command make all makes daplo.exe progsub.exe and progrw.exe The command make clean removes all the .o intermediate object files
Introducing a fitting library - we can try fitting a gaussian curveAgain we have the idea of a standard (or set of standard routines) to read in various types of data, and sets of calculation rouines to fit to the data using the library.
The MAIN program serves primarily to label the fitting parameters and select the reading and calculating routines.
PROGRAM FTEST C***** BASIC EXAMPLE ROUTINE FOR FITFUN C***** Link FTEST,MYREAD,GAUSS,FITFUN/LIB,LIBPG/LIB C HP setenv LPATH /lib:/usr/lib:/usr/local/lib:/usr/lib/X11R5 C HP f77 -o ftst ftst.o myread.o gauss.o libfitfun.a -lpgplot -lX11 -lU77 C IRIX f77 -o ftst ftst.o myread.o gauss.o libfitfun.a libpgplot.a -lX11 C MinGW f77 -o ftst ftst.o myread.o gauss.o libffn66m.a libpggw520.a \ C -mconsole -mwindows EXTERNAL MYREAD,GAUSS C C***** MYREAD IS A USER SUPPLIED ROUTINE TO READ IN ONE SPECTRUM EACH CALL C C***** GAUSS IS A USER'S ROUTINE TO CALCULATE THE REQUIRED FUNCTION C COMMON/TITLES/NAMES(40),TX,TY COMMON/TITLEP/NPARAS COMMON/WORK/W(30660) C***** THIS IS SUFFICIENT FOR 1600 DATA AND 10 PARAMETERS ! COMMON/VERSION/VERP C***** VERP IS TYPED BY THE PROGRAM INDICATING THE CURRENTLY RUNNING VERSION C CHARACTER*8 NAMES CHARACTER*4 PNAM CHARACTER*20 TX,TY C DATA NAMES/'CENTRE ','SIGMA','FLAT BAK','SLOPE', 1'HEIGHT',35*' '/ DATA PNAM/'ftst'/ DATA TX/'CHANNEL'/ DATA TY/'COUNTS'/ DATA NPARAS/5/ C VERP=2.0 CALL FITFUN(PNAM,MYREAD,GAUSS) C END The calculation of a simple gaussian: SUBROUTINE GAUSS(NPAR,PARM,NFIT,XUSE,YUSE,YRUSE,YCALC,F) C C CALC OF F AND YCALC FOR GAUSSIAN C DIMENSION PARM(NPAR) DIMENSION YCALC(NFIT),XUSE(NFIT),YUSE(NFIT) 1,YRUSE(NFIT),F(NFIT) C C The first parameter PARM(1) is the peak centre C PARM(2) is the peak deviation C PARM(3) is the flat background C PARM(4) is the sloping background C PARM(5) is the peak height C There are NFIT data values of Xuse Yuse and Yerror DO 10 I=1,NFIT C***** YCALC WILL CONTAIN CALCULATED FUNCTION FOR PLOTTING/LISTING YCALC(I)= PARM(5) * 1 EXP(-0.5*((XUSE(I)-PARM(1))/PARM(2))**2) 1 +PARM(3) + PARM(4)*XUSE(I) F(I)=YCALC(I)-YUSE(I) C***** F WILL BE MINIMISED BY the FITFUN LIBRARY by varying PARM(1-5) 10 CONTINUE C C RETURN END subroutine myread(in,xx,yy,yr,ntext) c****** routine to read from a simple GSAS data file c xx out first column value c yy out second column value c yr out third column value c in out number of lines read c inmax in maximum lines expected c R. Ghosh, Ill, October 2004 parameter (inmax=1600) real xx(inmax),yy(inmax),yr(inmax) common/infile/fname character*100 fname c***** the filename to be read has a maximum length of 100 characters character ntext*(*) write(6,1) 1 format(' Give input filename : ',$) read(5,2) fname 2 format(a) open(unit=2,file=fname) c***** read title read(2,5) ntext 5 format(a) read(2,6) read(2,6) c***** skip 2 lines 6 format(1x) do i=1,inmax,5 read(2,*,end=20) (yy(n),yr(n),n=i,i+4) j=i end do 20 continue close(unit=2) do i=1,j xx(i)=i*0.1-2. end do c***** not required further - close input write(6,3) j 3 format(i3,' lines have been read in') in=j return endWe can add another few lines to the makefileftest.exe: ftest.f myread.o gauss.o g77 -o ftest ftest.f myread.o gauss.o libffn66m.a libpggw520.a \ -wl,subsystem,console -mwindows myread.o: myread.f g77 -c myread.f gauss.0: gauss.f g77 -c gauss.fand then use the command gmake ftest.exe then test it after reading in the data file S100.gsa
A more complicated case allows one to fit up to four peaks using different model functions but sharing the read-in routineff4.exe: ff4.f calskp.f calskpx.f g77 -o ff4 ff4.f myread.o calskp.f calskpx.f libffn66m.a libpggw520.a \ -mconsole -mwindows
These programs test data and libraries have been collated in a zip file together with a C version of the PGPLOT library and its header file cpgplot.h.