C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C VIEWer for Miscellaneous FTMW spectra C ---- - C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c This program allows viewing of multiple spectra: c c I. From a .FAR spectral archive as produced by FILMAN c II. An assortment of individual spectral files c as specified by the file LIST.DAT or similar produced with programs c FFTLIST or FILMAN or FFT8 in AUTO-SCAN mode c III. Any spectral files found in the local directory c c c Interdependence of displays: c c SUMMAR-----LOOKSP c | c --------INTCOM------LOOKSP c c where: SUMMAR - summary screen of fringe voltages and interferogram c amplitudes as a function of frequency c INTCOM - interferogram comparisons, nine interferograms per screen c LOOKSP - individual interferogram and its FFT c c NOTE: these routines, though they have the same names, are c different and incompatible with those in program VIEW for c viewing AUTOscanned spectra c C Ver 3.I.2010 ---- Zbigniew KISIEL ---- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C c 4.05.08: Creation from V6 and VKIEL c 11.05.08: Synchronisation with FTMW_20 c 17.05.08: More synchronisation: modified file standard, flexible if, c synthetic spectra c 28.05.08: cleaning up c 5.11.08: subtraction of selected line from the interferogram c 14.01.09: interferogram points concatenated into a single vector + mods c 31.12.09: display comments and G option in SUMMAR, CTRL I in LOOKSP C 3.01.10: port to IVF C C_____________________________________________________________________________ c c MAXSPE - upper limit on the number of spectra that can be read in c MAXPTS - upper limit on the number of points in each interferogram c NIVOLS - number of diagnostic interferogram range voltages stored c in SCAN.DAT c MAXALL - upper limit on the total number of concatenated interferogram c points (dynamic) c c INTALL - vector containing data points of all interferograms c IPOINT(i) - pointer to interferogram i in order of frequency c INTPT(i,j) - indices of points in INTALL for interferogram i, c j=1 is the first point, j=2 is the last point c FREQ(i) - centre frequency for interferogram i, these are sorted c VSTEPS(i) - voltage spacing per pixel for interferogram i (Volts) c TSTEPS(i) - time spacing per pixel for interferogram i (seconds) c NREPS(i) - number of recorded points for interferogram i c NSKIP(i) - no of points to be skipped at the beginning for interferogram i c NSKIP1(i) - no of points to be skipped at the end for interferograms i c NAVE(i) - no averages for interferogram i c FNAMS(i) - file name for interferogram i c NDETS(i) - detection type index for interferogram i c FIKHZS(i) - intermediate frequency for interferogram i c ISEEN(i) - colour index used to indicate whether the interferogram has c already been looked at (11) or not (15) c COMNT(i) - comment string for interferogram i c TIMD(i) - time and date string for interferogram i c SAMPL(i) - sample string for interferogram i c DETVOL(i,j) - two cavity fringe voltages (j=1,2) for interferogram i written c after the data points when in autoscanning mode (indicated by c -ve value of naver) c c IDATA(i) - data points for interferogram i c FCENT - centre frequency for the current interferogram c VSTEP - voltage spacing per pixel for the current interferogram (Volts) c TSTEP - time spacing per pixel for the current interferogram (seconds) c NREP - number of recorded points for the current interferogram c (truncated to MAXPTS if longer) c NSKIPS - no of points to be skipped at the beginning for the current infgram c NSKIPE - no of points to be skipped at the end for the current interferogram c NAVER - no averages for the current interferogram c FILNAM - file name for the current interferogram c NDET - detection type for the current interferogram c FIFKHZ - intermediate frequency for the current interferogram c c VOLINT - interferogram amplitudes for various values of NSKIPS (the c central value is for NSKIPS as set during acquisition) C_____________________________________________________________________________ c c Colours: c C 0 - black 4 - red 8 - dark grey 12 - light red C 1 - blue 5 - magenta 9 - light blue 13 - light magenta C 2 - green 6 - brown 10 - light green 14 - yellow C 3 - cyan 7 - white 11 - light cyan 15 - bright white C C----------------------------------------------------------------------------- C C O M P I L A T I O N: C----------------------------------------------------------------------------- C C This version will only compile satisfactorily with C Compaq Visual Fortran 6.50 (and possibly with not too distant earlier C versions of Microsoft Powerstation Fortran) C Note that there is compatibility with older versions in that C C Compilation is now to be for QWIN graphics - this necessitates explicit C programming out of several unnecessary frills, but results in smoother C launch of the program than is possible with the STANDARD graphics as used C previously. C C-------------------------------- C Command line compilation: C-------------------------------- C C Simplest compilation for the local machine: C C df -static -libs=qwin -fpscomp:filesfromcmd viewm.for C C Optimised compilation for any PENTIUM: C C df -nodebug -traceback -arch=pn1 -tune=pn1 C -fast -static -libs=qwin -fpscomp:filesfromcmd viewm.for C C Other processor options are pn2,pn3,pn4,k6_2,k7 C C-------------------------------- C Visual Studio compilation: C-------------------------------- C C FORTRAN: /compile_only /fpscomp:filesfromcmd C /libs:qwin /nologo /nopdbfile /optimize:3 /traceback /tune:pn1 C /architecture:pn1 /static C C LINK: kernel32.lib /nologo /subsystem:windows /pdb:none C /machine:IX86 /out:"Debug/v6.exe" C C The use of /check:all FORTRAN option is also recommended, but only for C debugging. C C-------------------------------------------------------- C Visual Studio Net compilation for Intel Visual Fortran: C-------------------------------------------------------- C C FORTRAN: /nologo /O3 /G5 /fpscomp:filesfromcmd /warn:unused C /Qsave /module:"Debug/" /object:"Debug/" /traceback C /libs:qwin /c C C LINK /OUT:"Debug/v6_20.exe" /NOLOGO /SUBSYSTEM:WINDOWS C /MACHINE:IX86 kernel32.lib C C Note: the -Qsave compiler option is equivalent to the previous -static C compiler option for static variable allocation. The -static option C now has the meaning of preventing linking with shared libraries. C C----------------------------------------------------------------------------- C S T A R T U P: C----------------------------------------------------------------------------- C C Program startup: C C 1/ in Win95,Win98 just call the program from the command line set to C the directory containing the data C 2/ in Win2000 use the command: C START /d . c:\fft\viewm.exe C from the command line set to the directory containing the data, C assuming the PATH leads to the directory containing VIEWM.EXE C C///////////////////////////////////////////////////////////////////////////// C C C Declarations associated with storage of concatenated interferograms in C a single vector C MODULE globdat c c INTALL - vector with points from all interferogram c MAXALL - size of INTALL, expanded dynamically when necessary c IALL - counter of used points in INTALL c I4TEMP - temporary buffer used when increasing MAXALL c INTPT - indices of first and last points in INTALL for a given c interferogram c INDXI - index used for transferring points from INTALL to IDATA c INDXJ - as above c parameter (maxspe=550,maxpts=25000) INTEGER*4 intall,maxall,iall,i4temp,intpt(maxspe,2),indxi,indxj ALLOCATABLE :: intall(:),i4temp(:) c END MODULE globdat C C///////////////////////////////////////////////////////////////////////////// C C M A I N P R O G R A M c c C...Initialization commands for graphics. The three structured C variables contain coordinates: C curpos.row and curpos.col - cursor coordinates (INTEGER*2) C ixy.xcoord and ixy.ycoord - pixel coordinates (INTEGER*2) C wxy.wx and wxy.wy - window coordinates (REAL*8) C C USE DFLIB use globdat C RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy C_____________________________________________________________________________ c parameter (nivols=7,maxsmo=1001) PARAMETER (ntextc=0, ntextb=7) c character fnams(maxspe)*12,filnam*30,filarc*30,dirnam*50 real*4 detvol(maxspe,2),volint(maxspe,nivols),spol(maxsmo) integer*4 idata(maxpts),nave(maxspe), * ioldat(maxpts),itemp(maxpts),dummy4 integer*2 iseen(maxspe),ipoint(maxspe),dum2,ndets(maxspe) INTEGER*2 maxx,maxy,LINOFS,mymode,myrows,mycols real*8 fifmhz,fifkhz,fifmho,fifkho,fikhzs(maxspe) real*8 freq(maxspe) integer*2 nreps(maxspe),nskip(maxspe),nskip1(maxspe),dummy2 real*4 tsteps(maxspe),vsteps(maxspe) character comnt(maxspe)*50,sampl(maxspe)*20,timd(maxspe)*20 c common /scans/freq,wmult,ipoint,filarc common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /scans2/fikhzs,detvol,volint,iseen,nscans,ndets common /sfiles/fnams common /scan3/nreps,nskip,nskip1,tsteps,vsteps,nave COMMON /SCTEXT/COMNT,SAMPL,TIMD common /smooth/ioldat,itemp,spol COMMON /limits/wxy,maxx,maxy,LI NOFS,curpos,ixy, * mymode,myrows,mycols common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto c irefr=0 mxnrep=0 mnnrep=maxpts c call startg(iconf) <----- dummy4=passdirkeysqq(.true.) c C...HEADER C numfonts = INITIALIZEFONTS ( ) fontnum = SETFONT ('t''Arial''h75w25ei') dummy4=setbkcolor(ntextb) call clearscreen($gclearscreen) c NBOTL=120 dummy=setcolor(int2(15)) CALL MOVETO (INT2( 0), INT2(NBOTL), ixy) dummy=lineto(INT2( maxx), INT2(NBOTL)) CALL MOVETO (INT2( 0), INT2(0), ixy) dummy=lineto(INT2( maxx), INT2(0)) dummy=setcolor(int2(8)) dummy=floodfill(int2(1),int2(1),int2(15)) c nvert=(NBOTL-75)/2 nhor=28 c dummy=setcolor(int2(11)) CALL MOVETO (INT2(nhor), INT2(nvert), ixy) CALL OUTGTEXT('V6_20 -') dummy=setcolor(int2(9)) CALL MOVETO (INT2(nhor+1), INT2(nvert+1), ixy) CALL OUTGTEXT('V6_20 -') dummy=setcolor(int2(1)) CALL MOVETO (INT2(nhor+2), INT2(nvert+2), ixy) CALL OUTGTEXT('V6_20 -') c dummy=setcolor(int2(11)) fontnum = SETFONT ('t''Arial''h60w20ei') CALL MOVETO (INT2(nhor+220), INT2(nvert+11), ixy) CALL OUTGTEXT('Viewer for FTMW Spectra') dummy=setcolor(int2(9)) CALL MOVETO (INT2(nhor+220+1), INT2(nvert+11+1), ixy) CALL OUTGTEXT('Viewer for FTMW Spectra') dummy=setcolor(int2(1)) CALL MOVETO (INT2(nhor+220+2), INT2(nvert+11+2), ixy) CALL OUTGTEXT('Viewer for FTMW Spectra') c dummy=setcolor(int2(15)) CALL MOVETO (INT2( 0), INT2(NBOTL+ 32), ixy) dummy=lineto(INT2( maxx), INT2(NBOTL+ 32)) dummy=setcolor(int2(8)) dummy=floodfill(int2(1),INT2(NBOTL+30),int2(15)) c dummy=setcolor(int2(0)) CALL MOVETO (INT2( 0), INT2(NBOTL+ 32), ixy) dummy=lineto(INT2( maxx), INT2(NBOTL+ 32)) dummy=setcolor(int2(7)) CALL MOVETO (INT2( 0), INT2(NBOTL+ 1), ixy) dummy=lineto(INT2( maxx), INT2(NBOTL+ 1)) c fontnum = SETFONT ('t''Arial''h20w10') dummy=setcolor(int2( 0)) CALL MOVETO (INT2( 11), INT2(NBOTL+ 7), ixy) CALL OUTGTEXT('version 3.I.2010') CALL MOVETO (INT2( maxx-169), INT2(NBOTL+ 7), ixy) CALL OUTGTEXT('Zbigniew KISIEL') dummy=setcolor(int2(15)) CALL MOVETO (INT2( 10), INT2(NBOTL+ 6), ixy) CALL OUTGTEXT('version 3.I.2010') CALL MOVETO (INT2( maxx-170), INT2(NBOTL+ 6), ixy) CALL OUTGTEXT('Zbigniew KISIEL') c nrlin=nint(real(NBOTL+32)/(real(maxy)/real(myrows)))+2 call settextposition(nrlin,44,curpos) dummy=settextcolor(int2(ntextc)) write(*,'(''up to'',i6,'' spectra'',$)')maxspe call settextposition(nrlin+1,44,curpos) write(*,'(''up to'',i6,'' points in each spectrum'')')maxpts c C...Warn of missing config file C 70 call settextposition(11,1,curpos) if(iconf.eq.0)then dummy=setcolor(int2(12)) fontnum = SETFONT ('t''Arial''h18w9e') CALL MOVETO (INT2(100), INT2(220), ixy) CALL OUTGTEXT( * 'Configuration file C:\FFT\V6.CFG was not found:') CALL MOVETO (INT2(250), INT2(240), ixy) fontnum = SETFONT ('t''Arial''h18w9i') CALL OUTGTEXT( * 'default sized window of 800x540 pixels will be used') call settextposition(15,1,curpos) endif c c c...Decide on type of input c maxall=100000 ALLOCATE(intall(maxall),stat=ialloc) c 555 call inpout(iarch,filarc,iexit,dirnam,ntsys) <----- if(iexit.eq.1)then dum2=setexitqq(qwin$exitnopersist) stop endif c c- - - - - - IARCH=0 - read individual spectral files - - - - - - - - - - - - c c either as found in current directory or as listed in c a listing file (from FILMAN, SVIEW, FFTLIST) c c...read LIST.DAT until EOF, errors in input skip the line in question so that c comment lines are allowed, but best be initiated with some nonnumeric c symbol such as $ or ! c 557 if(iarch.eq.0)then c if(irefr.eq.1)goto 558 c open(2,file='list.dat',status='old',err=1101) write(*,1102)'Using the standard listing file LIST.DAT' 1102 format(1x/' ----> ',a) write(filarc,'(a)')'Spectra listed in LIST.DAT' goto 102 1101 write(*,1102)'Usable LIST.DAT not found' c 101 write(*,100) 100 format(1x/' OPTIONS:'/ * ' ENTER = read spectral files in current directory'/ * ' file name = name of listing file for input of spectra'/ * ' minus sign = exit'//25x,'..... ',$) read(*,'(a)',err=101)filnam c c...exit c if(filnam(1:1).eq.'-')then dummy=setexitqq(qwin$exitnopersist) stop endif c c - - - - - - read spectra as found in the current directory c 558 if(filnam(1:1).eq.' '.or.filnam(1:1).eq.char(0))then iarcht=-1 call inpspe(nfil,ntsys) <----- if(nfil.eq.0)then write(*,1115) 1115 format(1x//' ----> The directory ',$) DUMmy2=settextcolor(int2(12)) n=len_trim(dirnam) write(*,'(1x,a)')dirnam(1:n) DUMmy2=settextcolor(int2(ntextc)) write(*,'(8x, * ''appears to contain NO FTMW spectral files'')') goto 101 endif c mxnrep=0 nscans=1 do 1006 n=1,nfil filnam=fnams(n) call readsp(iarch,nscans,filnam,iread) <----- c if(iread.gt.0)then ISEEN(nscans)=15 nreps(nscans)=nrep nskip(nscans)=nskips nskip1(nscans)=nskipe tsteps(nscans)=tstep vsteps(nscans)=vstep nave(nscans)=naver fnams(nscans)=fnams(n) if(nrep.gt.mxnrep)mxnrep=nrep ipoint(nscans)=NSCANS ndets(nscans)=ndet fikhzs(nscans)=fifkhz c indxi=intpt(nscans,1)-1 do 1010 i=1,nrep indxi=indxi+1 idata(i)=intall(indxi) ioldat(i)=idata(i) 1010 continue c call baksub(51) <----- c do 1011 i=1,nivols mindat= 1400000000 maxdat=-1400000000 j=nskips-(4-i)*20 do 1012 jj=j,nrep jjj=jj if(jjj.lt.1)jjj=1 if(idata(jjj).lt.mindat)mindat=idata(jjj) if(idata(jjj).gt.maxdat)maxdat=idata(jjj) 1012 continue volint(nscans,i)=1000.*vstep*(real(maxdat)-real(mindat)) 1011 continue c nscans=nscans+1 endif 1006 continue c nscans=nscans-1 nrep=mxnrep write(filarc,'(a)')'Spectra in local directory' c if(nscans.eq.0)then write(*,1115) DUMmy2=settextcolor(int2(12)) n=len_trim(dirnam) write(*,'(1x,a)')dirnam(1:n) DUMmy2=settextcolor(int2(ntextc)) write(*,'(8x, * ''appears to contain NO FTMW spectral files'')') goto 101 endif c c write(*,'(1x//i5,'' spectra have been identified''/)')nscans c write(*,*)(fnams(j),j=1,nscans) c pause c goto 1105 c endif c c - - - - - - use the specified listing file c open(2,file=filnam,status='old',err=101) filarc='Spectra listed in '//filnam iarcht=0 c 102 nscans=0 c 1 if(nscans.ge.maxspe)then read(2,3,end=2,err=1)filnam,vstep, * nrep,tstep,vstep,(vstep,j=1,nivols) if(vstep.eq.0.)goto 1 nscans=nscans+1 goto 1 else read(2,3,end=2,err=1)fnams(nscans+1),freq(nscans+1), * nrep,tstep,vstep,(volint(nscans+1,j),j=1,nivols) if(volint(nscans+1,nivols).eq.0.)goto 1 endif 3 format(a12,f12.3,i6,2f6.3,7f10.3) nscans=nscans+1 ipoint(nscans)=NSCANS if(nrep.gt.mxnrep)mxnrep=nrep if(nrep.lt.mnnrep)mnnrep=nrep goto 1 c 2 close(2) c write(*,4)NSCANS,mnnrep,mxnrep 4 FORMAT(1x//' The definition file has',i6,' interferograms'/ * ' minimum length = ',i6/ * ' maximum length = ',i6/) if(mxnrep.gt.maxpts)then dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) write(*,3470)maxpts dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) endif if(nscans.eq.0.or.mnnrep.lt.50.or.mxnrep.lt.50)then write(*,1102)'List file rejected, try again' goto 101 endif if(nscans.gt.maxspe)then dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) nscans=maxspe write(*,3471)nscans dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) endif 3470 format(' ***** Interferograms will be chopped to ', * i6,' points') 3471 format(' ***** Only the first ',i5, * ' interferograms will be read') write(*,3473) 3473 format(1x/50x,'Press E N T E R ',$) read(*,'(i1)',err=3472)j 3472 call clearscreen($gclearscreen) write(*,'(1x/'' N O W R E A D I N G:''/)') c mxnrep=0 do 5 n=1,nscans filnam=fnams(n) call readsp(iarch,n,filnam,iread) <----- ISEEN(N)=15 nreps(n)=nrep nskip(n)=nskips nskip1(n)=nskipe tsteps(n)=tstep vsteps(n)=vstep nave(n)=naver if(nrep.gt.mxnrep)mxnrep=nrep ndets(n)=ndet fikhzs(n)=fifkhz 5 continue nrep=mxnrep c c...SORT interferograms in frequency c 1105 if(nscans.gt.1)CALL SORTH <----- c endif c c- - - - - - IARCH=1 - read files from spectral archive - - - - - - - - - - - c 556 if(iarch.eq.1)then iarcht=1 open(3,file=filarc,form='binary',status='old') c call clearscreen($gclearscreen) write(*,'(1x/'' N O W R E A D I N G:''/)') n=0 c c...main loop for extraction of spectra from archive c ichop=0 411 call readsp(iarch,n,filnam,iread) <----- if(iread.eq.-1)goto 410 if(iread.eq.2)ichop=1 ipoint(n)=n fnams(n)=filnam nreps(n)=nrep nskip(n)=nskips nskip1(n)=nskipe tsteps(n)=tstep vsteps(n)=vstep nave(n)=naver if(nrep.gt.mxnrep)mxnrep=nrep if(nrep.lt.mnnrep)mnnrep=nrep ndets(n)=ndet fikhzs(n)=fifkhz c c...subtract baseline using a 51 point smooth indxi=intpt(n,1)-1 do 210 i=1,nrep indxi=indxi+1 idata(i)=intall(indxi) ioldat(i)=idata(i) 210 continue call baksub(51) <----- c c...determine diagnostic voltages do 201 i=1,7 mindat= 1400000000 maxdat=-1400000000 j=nskips-(4-i)*20 do 200 jj=j,nrep jjj=jj if(jjj.lt.1)jjj=1 if(idata(jjj).lt.mindat)mindat=idata(jjj) if(idata(jjj).gt.maxdat)maxdat=idata(jjj) 200 continue volint(n,i)=1000.*vstep*(real(maxdat)-real(mindat)) 201 continue c if(n.eq.maxspe)then dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) write(*,412)maxspe dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) 412 format(1x/' ***** The maximum of',i5,' interferograms', * ' reached - no more will be read',$) else goto 411 endif 410 close(3) if(ichop.eq.1)then dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) write(*,454)maxpts dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) endif 454 format(1x/' ***** Interferograms have been chopped to',i6, * ' points',$) nscans=n nrep=mxnrep c write(*,540)NSCANS,filarc,mnnrep,mxnrep,iall 540 FORMAT(1x//i6,' interferograms have been read from ',a/ * ' minimum length = ',i6/ * ' maximum length = ',i6/ * ' total points =' ,i7/) write(*,3473) read(*,'(i1)',err=3475)j c 3475 if(nscans.gt.1)CALL SORTH <----- c do 203 n=1,nscans ISEEN(N)=15 203 continue c endif c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...plot statistics c iarch=iabs(iarcht) call summar(iarch,irefr,dirnam) <----- if(irefr.eq.1)then if(filarc(1:16).eq.'Spectra in local')then iarch=0 filarc=' ' endif if(iarch.eq.1)then goto 556 else filnam=' ' goto 557 endif endif c stop end c c--------------------------------------------------------------------------- c subroutine readsp(iarch,nspe,filnam,iread) c c...Read spectrum recorded in the old style and also in version with FIFKHZ c introduced 12.V.2008 c c IARCH=0 Read spectral file FILNAM and store it as spectrum number NSPE c - if there are more than MAXPTS data points they are truncated c to MAXPTS c IARCH=1 Read spectrum from the currently open spectral archive, c reading past the end of the archive results in IREAD=-1 c c IREAD=-1 end of file reached while attempting to read spectrum c IREAD= 0 spectrum could not be read properly c IREAD= 1 spectrum read in without problems c IREAD= 2 spectrum chopped to maxpts c USE DFLIB use globdat PARAMETER (nivols=7) PARAMETER (ntextc=0, ntextb=7) c real*8 fifmhz,fifkhz,fifmho,fifkho,fikhzs(maxspe) real*8 fcent,freq(maxspe) character timdat*20,coment*50,sample*20,INTEXT*30 character filnam*30,cdummy*6,filarc*30 character comnt(maxspe)*50,sampl(maxspe)*20,timd(maxspe)*20 integer*4 idata(maxpts),idummy integer*2 ipoint(maxspe),ndets(maxspe) real*4 detvol(maxspe,2),volint(maxspe,nivols) c common /scans/freq,wmult,ipoint,filarc common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /scans2/fikhzs,detvol,volint,iseen,nscans,ndets COMMON /SCTEXT/COMNT,SAMPL,TIMD common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto c iread=1 if(iarch.eq.0)then write(*,'(4x,a12,$)')filnam(1:12) OPEN(3,FILE=FILNAM,FORM='BINARY',ERR=503,STATUS='OLD') else read(3,end=510)cdummy,filnam(1:12),cdummy write(*,'(4x,a12,$)')filnam(1:12) endif C c Data file is binary but all header information is ASCII so that it can c be inspected with an editor, and modified with a hex editor (THE NUMBERS C OF CHARACTERS ARE GIVEN IN BRACKETS) c c SAMPLE (20) = sample name c NREP (5) = number of points in the interferogram c NSKIPS (5) = number of points from beginning of interf. to be skipped c for FFT c NSKIPE (5) = number of points from end of interf. to be skipped c for FFT c FCENT (12) = excitation frequency (MHz) c FIFKHZ (10) = intermediate frequency (kHz) which is a signed quantity c in F_sweeper=F_synt+F_if c NDET (3) = detection type: c 0 = signal demixed to zero IF in second mixer c 1 = no demixing in second mixer, signal is on IF c TSTEP (10) = horizontal spacing between points (seconds) c VSTEP (10) = vertical spacing between points (volts) c NAVER (7) = number of averages c COMENT (50) = comment c TIMDAT (20) = time and date c c This block is followed by NREP data points from IDATA, 4-bytes each. c If NAVER is negative (AUTO-Scan) then VFRPK and VFRRET follow IDATA, c also 4-byte long each c c NOTE: FIFKHZ and NDET are have been introduced on 12.V.2008. Previous c spectra are for FIFKHZ=20000. and NDET=0, and are recognised by c the small value of TSTEP, which is seen on input in place of c FIFKHZ. In that case intermediate frequency is set to zero c and the order of parameters after FCENT is: C c TSTEP (10) = horizontal spacing between points (seconds) c VSTEP (10) = vertical spacing between points (volts) c NAVER (7) = number of averages c COMENT (50) = comment c TIMDAT (20) = time and date C c READ(3,end=510)sample READ(3,end=510)intext(1:27) READ(intext,'(3i5,f12.5)',err=510)nrep,nskips,nskipe,fcent if(nrep.lt.50)goto 510 if(fcent.lt.500.or.fcent.gt.26000.)goto 510 c c...read the 10 byte ASCII number immediately after FCENT and decide whether c this is c 1/ TSTEP = file in old standard c 2/ FIFKHZ = file in new standard c READ(3)intext(1:10) READ(intext(1:10),*)tstep if(abs(tstep).gt.0.1d0)then fifkhz=tstep fifmhz=fifkhz*0.001D0 read(3)intext(1:13) read(intext,'(i3,1pe10.3)')ndet,tstep else fifkhz=0.d0 fifmhz=0.d0 ndet=0 endif READ(3)intext(1:17) c READ(intext,'(1pe10.3,i7)')vstep,naver READ(3,end=510)coment,timdat c IF(IARCH.EQ.1)NSPE=NSPE+1 freq(nspe)=fcent COMNT(NSPE)=coment timd(nspe)=timdat sampl(nspe)=sample c c...Expand the buffer for interferograms if necessary c if(iall+nrep.gt.maxall)then ALLOCATE(i4temp(maxall),stat=ialloc) do 21 n=1,iall i4temp(n)=intall(n) 21 continue DEALLOCATE(intall,stat=ialloc) maxall=maxall*2 ALLOCATE(intall(maxall),stat=ialloc) do 22 n=1,iall intall(n)=i4temp(n) 22 continue DEALLOCATE(i4temp,stat=ialloc) endif C c...Read interferogram data points c if(nspe.eq.1)iall=0 DO 20 N=1,Nrep READ(3,end=510)idummy if(n.le.maxpts)then iall=iall+1 intall(iall)=idummy endif 20 CONTINUE c if(nrep.gt.maxpts)then iread=2 nskipe=nskipe-(nrep-maxpts) if(nskipe.lt.1)nskipe=1 nrep=maxpts if(nrep-nskips-nskipe.le.5)then nskipe=5 nskips=5 endif endif if(nspe.eq.1)then intpt(nspe,1)=1 intpt(nspe,2)=nrep else intpt(nspe,1)=intpt(nspe-1,2)+1 intpt(nspe,2)=intpt(nspe-1,2)+nrep if(intpt(nspe,2).ne.iall)then write(*,'(1x//'' PROBLEM: iall='',i8,'' inpt(nspe,2)='', * i8/)')iall,intpt(nspe,2) pause stop endif endif c if(naver.lt.0)then read(3,end=510)detvol(nspe,1) read(3,end=510)detvol(nspe,2) else detvol(nspe,1)=0.0 detvol(nspe,2)=0.0 endif C if(iarch.eq.0)CLOSE(3) goto 501 c 503 dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) write(*,504)filnam 504 format(1x/' ***** Cannot open file: ',a,$) dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) write(*,'(1x)') iread=0 return 510 iread=-1 c 501 return end c c--------------------------------------------------------------------------- c subroutine summar(iarch,irefr,dirnam) c C Routine to display summary information on all acquired spectra and allow c cursor selection of spectrum/spectra for display c c IARCH - when 0 then input has been from individual data files c when 1 then input has been from an archive c IREFR - refresh flag, set by SUMMAR. Normally zero but when set to 1 c on exit signals M/P to reread data files c DIRNAM - character variable with the name of the directory c c containing files c c VRANGE = voltage range (mV) for scaling of interferogram ranges c DRANGE = voltage range (mV) for scaling of fringe voltages c c...declarations necessary for graphics c USE DFLIB use globdat c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy INTEGER*2 maxx,maxy,LINOFS,dummy,inkey,mymode,myrows,mycols INTEGER*4 dummy4 logical*2 true real*8 wmin,wmax,fstart,fend,fmark,rint,wrange,fincr,freqv,fp REAL*8 YYSTEP,YSHIFT,RRSMAL,fchang real*8 pi,phase,rifph,rifamp,rifsh,riffr character kk,emplin*80,lwork1*80,outstr*21,dirnam*50 COMMON /limits/wxy,maxx,maxy,LINOFS,curpos,ixy, * mymode,myrows,mycols COMMON /plotda/wmin,wmax,RRSMAL,YYSTEP,YSHIFT COMMON /lines/emplin common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto c c...declarations for spectra c parameter (nivols=7, vrange=50.0,true=.true.) PARAMETER (ntextc=0, ntextb=7, nbordc=15, ncursc=14) PARAMETER (Nmaxpt=131072,maxsyn=1000000,maxpar=10) PARAMETER (pi=3.141592653589793d0) c character fnams(maxspe)*12,filarc*30,fspec*12 real*4 detvol(maxspe,2),volint(maxspe,nivols),p(nmaxpt) real*8 fifmhz,fifkhz,fifmho,fifkho,fikhzs(maxspe) integer*4 idata(maxpts),nave(maxspe) integer*2 iseen(maxspe),ipoint(maxspe),ndets(maxspe) real*8 freq(maxspe),fstep,df,bw,fs_st,fs_en,synspe(maxsyn),fstepm, * x0,x1,x2,cl0,cl1,cl2 integer*2 nreps(maxspe),nskip(maxspe),nskip1(maxspe) real*4 tsteps(maxspe),vsteps(maxspe) c real*8 xpfit(maxpts),ypfit(maxpts),sigpts(maxpts),acon(maxpar), * COVAR(maxpar,maxpar),ALPHA(maxpar,maxpar),CHISQ, * ALAMDA,ochisq integer*4 lista(maxpar) character comnt(maxspe)*50,sampl(maxspe)*20,timd(maxspe)*20 c EXTERNAL SUBIF,SUBLIN c common /scans/freq,wmult,ipoint,filarc common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /scans2/fikhzs,detvol,volint,iseen,nscans,ndets common /scan3/nreps,nskip,nskip1,tsteps,vsteps,nave common /sfiles/fnams common /specf/p,fstep,npts,NFFT,NCALL COMMON /SCTEXT/COMNT,SAMPL,TIMD c if(idispl.ne.-1)idispl=-1 c nmark=nscans/2+1 fmark=freq(nmark) c iptick=1 ifring=0 volt=0.0 n=nivols/2+1 drange=0.0 DO 201 I=1,nscans if(volt.lt.volint(i,n))VOLT=VOLINT(I,N) if(-detvol(i,1).gt.drange)drange=-detvol(i,1) if(-detvol(i,2).gt.drange)drange=-detvol(i,2) 201 CONTINUE if(drange.ne.0.0)ifring=1 wmult=volt/(vrange*0.9) wmulti=wmult c WRITE(emplin,'(80(1H ))') C c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c G R A P H I C S C C...Start up the graphics C itxt=(mycols-80)/2 C C...definition of graphics viewport: first pixel coordinates of viewport c then real coordinates for scales c 179 dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) call clearscreen($GCLEARSCREEN) call setviewport(int2(2),int2(2*LINOFS-2), * int2(maxx-2),int2(maxy-2*LINOFS+2)) <-vvvv c if(irefr.eq.1)then irefr=0 fstart=0.d0 fend=0.d0 endif if(fstart.eq.0.d0.and.fend.eq.0.d0)then fstart=freq(1) fend=freq(nscans) fstart=fstart-0.02*(fend-fstart) fend=fend+0.02*(fend-fstart) if(fend.eq.fstart)then fstart=fstart-2. fend=fend+2. endif fincr=1.04*(fend-fstart)/(maxx-4) endif c c...complete screen refresh takes place from here (note that unlike in c LOOKSP the limits FSTART and FEND are in MHz c 178 wrange=wmult*vrange wmin= 0.0d0 wmax= 1.05*WRANGE YYSTEP=(wmax-wmin)/500.d0 RRSMAL=wmin-13.d0*YYSTEP dummy=setwindow(TRUE,fstart,RRSMAL,fend,wmax) <-wwww dummy4=setbkcolor(1) DUMMY2=SETCOLOR(int2( nbordc )) call clearscreen($GVIEWPORT) DUMMY2=SETCOLOR(int2( 0 )) CALL moveto_w(fstart,wmin,wxy) dummy=lineto_w(fstart,wmax) dummy=lineto_w(fend,wmax) DUMMY2=SETCOLOR(int2( nbordc )) dummy=lineto_w(fend,wmin) dummy=lineto_w(fstart,wmin) c c...Plot range voltages using one of the two possible display schemes c c 1/ spectral scheme (join dots as function of frequency) C 2/ histogram scheme - i.e. a stick diagram c if(idispl.eq.1)then DUMMY=SETCOLOR(int2(12)) red do 102 n=1,nivols VOLT=VOLINT(1,N) IF(VOLT.GT.WRANGE)VOLT=WRANGE CALL moveto_w(freq(1),DBLE(volT),wxy) if(n.eq.nivols/2+1)dummy=setcolor(int2(15)) white if(n.gt.nivols/2+1)dummy=setcolor(int2(14)) yellow DO 103 I=1,nscans freqv=freq(I) if(freqv.lt.fstart.or.freqv.gt.fend)goto 103 VOLT=VOLINT(I,N) IF(VOLT.GT.wmax)VOLT=wmax dummy=lineto_w(freqv,dble(volT)) 103 CONTINUE 102 continue else do 100 I=1,nscans freqv=freq(I) if(freqv.lt.fstart.or.freqv.gt.fend)goto 100 c VOLT=VOLINT(I,1) IF(VOLT.GT.wmax)VOLT=wmax DUMMY=SETCOLOR(int2(12)) red CALL moveto_w(freqv,0.d0,wxy) dummy=lineto_w(freqv,dble(volT)) c VOLT=VOLINT(I,4) IF(VOLT.GT.wmax)VOLT=wmax DUMMY=SETCOLOR(int2(15)) white CALL moveto_w(freqv,0.d0,wxy) dummy=lineto_w(freqv,dble(volT)) c DUMMY=SETCOLOR(int2(14)) yellow VOLT=VOLINT(I,5) IF(VOLT.GT.wmax)VOLT=wmax CALL moveto_w(freqv,0.d0,wxy) dummy=lineto_w(freqv,dble(volT)) c VOLT=VOLINT(I,6) IF(VOLT.GT.wmax)VOLT=wmax CALL moveto_w(freqv-fincr,0.d0,wxy) dummy=lineto_w(freqv-fincr,dble(volT)) CALL moveto_w(freqv+fincr,0.d0,wxy) dummy=lineto_w(freqv+fincr,dble(volT)) c VOLT=VOLINT(I,7) IF(VOLT.GT.wmax)VOLT=wmax CALL moveto_w(freqv-2.d0*fincr,0.d0,wxy) dummy=lineto_w(freqv-2.d0*fincr,dble(volT)) CALL moveto_w(freqv+2.d0*fincr,0.d0,wxy) dummy=lineto_w(freqv+2.d0*fincr,dble(volT)) c 100 continue c endif c c...plot fringe voltages, if available c if(ifring.eq.1)then detscl=WRANGE/drange DUMMY=SETCOLOR(int2(3)) cyan do 105 n=1,2 j=0 DO 106 I=1,nscans freqv=freq(I) if(freqv.lt.fstart.or.freqv.gt.fend)goto 106 j=j+1 rint=-detvol(i,n)*detscl+0.1*detscl if(rint.gt.wrange)RINT=WRANGE if(j.eq.1)then CALL moveto_w (freqv,rint,wxy) else dummy=lineto_w(freqv,rint) endif 106 CONTINUE 105 continue endif dummy=setcolor(int2(15)) C C...marker scale C yshift=wmin call marsca(fstart,fend,0.0d0,0.0d0,-1.d0) <----- c c...cursor c dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) dummy=setwritemode($GPSET) CALL setlinestyle(int2(#FFFF)) c c...information lines c dummy4=setbkcolor(ntextb) dummy=settextcolor(int2(12)) call settextposition(myrows,int2(itxt+36),curpos) call outtext(filarc(1:30)) dummy=settextcolor(int2(ntextc)) WRITE(OUTSTR,'(A,F8.2)')'Yrange /mV:',WRANGE CALL settextposition(myrows,int2(itxt+2),curpos) CALL outtext(OUTSTR(1:19)) CALL settextposition(myrows,int2(itxt+70),curpos) CALL outtext('H = help') c 771 write(lwork1,'('' file: '',12x,'' frequency:'', * f12.4,'' MHz'')'),freq(nmark) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork1) dummy=settextcolor(int2(12)) CALL settextposition(1,int2(itxt+8),curpos) CALL outtext(fnams(nmark)) c dummy=settextcolor(int2(2)) CALL settextposition(2,24,curpos) call outtext(comnt(ipoint(nmark))) dummy=settextcolor(int2(ntextc)) c c...options loop: c A,S - shift screen window left/right c Q,E - zoom in/out in frequency c W,Z - increase/decrease vertical scale c K,L - move cursor left/right c - cursor beginning/end C c I - go to screen comparing interferograms C O - go to spectrum under the cursor C G - go to a specified spectrum c c U - generate spectrum from loaded interferograms c P - toggle display style c R - restore initial scaling c - exit to M/P to reread data c - exit from program c 77 IK=INKEY(N) KK=CHAR(IK) c c...exit SUMMAR to M/P to reread data files (with ESC) c if(iK.eq.27)then irefr=1 return endif c c...exit to individual interferogram and its FFT (with O) c 303 if(KK.eq.'O'.or.kk.eq.'o')then nmark1=ipoint(nmark) nrep=nreps(nmark1) nskips=nskip(nmark1) nskipe=nskip1(nmark1) vstep=vsteps(nmark1) tstep=tsteps(nmark1) iseen(nmark)=11 ndet=ndets(nmark1) fifkhz=fikhzs(nmark1) fifmhz=0.001d0*fifkhz call looksp(nmark,nmark1) <----- GOTO 179 endif c c...go to comparison of interferograms (with I) c if(KK.eq.'I'.or.kk.eq.'i')then call intcom(nmark,iarch,irefr) <----- if(irefr.eq.1)return fmark=freq(nmark) if(fmark.lt.fstart.or.fmark.gt.fend)then frange=fend-fstart fstart=fmark-0.5d0*frange fend=fmark+0.5d0*frange endif GOTO 179 endif c c...zoom-in in frequency (with E) c if(KK.eq.'E'.or.KK.eq.'e')then FRange=Fend-Fstart Fchang=0.25D0*FRange IF(KK.EQ.'e')Fchang=0.45d0*FRange Fstart=Fmark-Fchang Fend=Fmark+Fchang if(fend-fstart.lt.0.1D0)then fstart=fmark-0.05d0 fend=fmark+0.05d0 endif c if(fstart.lt.freq(1))fstart=freq(1)-0.1d0 if(fend.gt.freq(nscans))fend=freq(nscans)+0.1d0 c fincr=1.04*(fend-fstart)/maxx goto 178 endif c c...zoom-out in frequency (with Q) c if(KK.eq.'Q'.or.KK.eq.'q')then FRange=Fend-Fstart Fchang=0.5D0*FRange IF(KK.EQ.'q')Fchang=0.1d0*FRange Fstart=Fstart-Fchang Fend=Fend+fchang c if(fstart.lt.freq(1))fstart=freq(1)-0.1d0 if(fend.gt.freq(nscans))fend=freq(nscans)+0.1d0 c fincr=1.04*(fend-fstart)/maxx goto 178 endif C C...Shift screen window to left (with A) C if(KK.eq.'A'.or.KK.eq.'a')then FRange=fend-fstart Fchang=FRange*0.5D0 IF(KK.EQ.'a')Fchang=FRange*0.1D0 fstart=fstart-Fchang fend= fend -Fchang IF(fstart.LT.freq(1))THEN fstart=freq(1)-0.1d0*frange fend=fstart+frange ENDIF 150 if(fmark.gt.fend)then nmark=nmark-1 fmark=freq(nmark) goto 150 endif if(fmark.lt.fstart)then fstart=fmark-frange*0.5d0 fend= fmark+frange*0.5d0 endif fincr=1.04*(fend-fstart)/maxx goto 178 endif C C...Shift screen window to right (with S) C if(KK.eq.'S'.or.KK.eq.'s')then FRange=fend-fstart Fchang=FRange*0.5D0 IF(KK.EQ.'s')Fchang=FRange*0.1D0 fstart=fstart+Fchang fend= fend +Fchang IF(fend.gt.freq(nscans))THEN fend=freq(nscans)+0.1d0*frange fstart=fend-frange ENDIF 151 if(fmark.lt.fstart)then nmark=nmark+1 fmark=freq(nmark) goto 151 endif if(fmark.gt.fend)then fstart=fmark-frange*0.5d0 fend= fmark+frange*0.5d0 endif fincr=1.04*(fend-fstart)/maxx goto 178 endif c c...cursor left, (with K) c if(KK.eq.'K'.or.KK.eq.'k'.or.ik.eq.-75)then dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) nmark=nmark-1 if(KK.eq.'K')nmark=nmark-9 if(nmark.lt.1)nmark=1 130 fmark=freq(nmark) if(fmark.lt.fstart)then nmark=nmark+1 goto 130 endif CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) CALL setlinestyle(int2(#FFFF)) dummy=setwritemode($GPSET) goto 771 endif c c...cursor right, (with L) c if(KK.eq.'L'.or.KK.eq.'l'.or.ik.eq.-77)then dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) nmark=nmark+1 if(KK.eq.'L')nmark=nmark+9 if(nmark.gt.nscans)nmark=nscans 131 fmark=freq(nmark) if(fmark.gt.fend)then nmark=nmark-1 goto 131 endif CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) CALL setlinestyle(int2(#FFFF)) dummy=setwritemode($GPSET) goto 771 endif c c...marker beginning,end c if(ik.eq.-71)then <- home dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) nmark=1 140 fmark=freq(nmark) if(fmark.lt.fstart)then nmark=nmark+1 goto 140 endif CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) CALL setlinestyle(int2(#FFFF)) dummy=setwritemode($GPSET) goto 771 endif c if(ik.eq.-79)then <- end dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) nmark=nscans 141 fmark=freq(nmark) if(fmark.gt.fend)then nmark=nmark-1 goto 141 endif CALL moveto_w(FMARK,wmax,wxy) dummy2=lineto_w(FMARK,wmin) CALL setlinestyle(int2(#FFFF)) dummy=setwritemode($GPSET) goto 771 endif C C...restore original spectrum and scaling (with R) C IF(KK.EQ.'R'.OR.KK.EQ.'r')THEN wmult=wmulti nmark=nscans/2+1 fmark=freq(nmark) fstart=freq(1) fend=freq(nscans) fstart=fstart-0.02*(fend-fstart) fend=fend+0.02*(fend-fstart) if(fend.eq.fstart)then fstart=fstart-2. fend=fend+2. endif fincr=1.04*(fend-fstart)/(maxx-4) goto 178 ENDIF c c...change of vertical scaling - zoom out (with Z) c if(KK.eq.'Z'.or.kk.eq.'z')then wmult=wmult*1.1 if(kk.eq.'Z')wmult=wmult*1.364 goto 178 endif c c...change of vertical scaling - zoom in (with W) c if(KK.eq.'W'.or.kk.eq.'w')then wmult=wmult/1.1 if(kk.eq.'W')wmult=wmult/1.364 goto 178 endif C C...toggle display style - between spectrum and histogram C if(KK.eq.'P'.or.kk.eq.'p')then idispl=-idispl goto 178 endif c c...Help screen (with H) c if(KK.eq.'H'.or.KK.eq.'h')then dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) CALL clearscreen($GCLEARSCREEN) WRITE(*,232) 232 FORMAT(1x/'_____S U M M A R Y screen commands', * 45(1H_)/// * 9x,' W/Z - zoom-in/zoom-out in intensity'/ * 9x,' E/Q - zoom-in/zoom-out in frequency'/ * 9x,' K/L - cursor left/right'/ * 9x,' A/S - window left/right'/ * 9x,' - cursor to beginning/end'// * 9x,' caps on/off - fast/slow change in the above'// * 9x,' U - generate synthetic spectrum'/ * 9x,' P - toggle spectrum <-> histogram display style'/ * 9x,' R - restore initial settings'// * 9x,' I - go to comparison of interferograms'/ * 9x,' O - go to spectrum under the cursor'/ * 9x,' G - go to a specified spectrum'// * 9x,' - (followed by Y) exit from the program'/ * 9x,' - reread input data'/// * 9x,' Press ENTER to exit this HELP screen'/) 107 IK=INKEY(J) IF(IK.NE.13)GOTO 107 CALL clearscreen($GCLEARSCREEN) CALL settextposition(2,int2(itxt+1),curpos) call outtext(emplin) GOTO 178 endif c c...Go to a specified spectrum (with G) c if(KK.eq.'G'.or.KK.eq.'g')then dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) dummy=displaycursor($gcursoron) 300 CALL clearscreen($GCLEARSCREEN) c nsypt=1.d0+(fs_en-fs_st)/df WRITE(*,301) 301 FORMAT(1x/'_____G o t o a s p e c i f i e d ', * ' s p e c t r u m',23(1H_)//) WRITE(*,304) 304 FORMAT(1x/ * ' File name (ENTER to exit) = ',$) c read(*,'(a)',err=300)fspec c do 302 i=1,nscans if(fnams(i).eq.fspec)then KK='O' nmark=i fmark=freq(nmark) frange=fend-fstart fstart=fmark-0.5d0*frange fend=fmark+0.5d0*frange if(fstart.lt.0.d0)then fstart=0.d0 fend=fstart+frange endif dummy=displaycursor($gcursoroff) goto 303 endif 302 continue c C...file not found, return to main display c call BEEPQQ(2000,250) dummy=displaycursor($gcursoroff) CALL clearscreen($GCLEARSCREEN) CALL settextposition(8,int2(itxt+1),curpos) call outtext(emplin) GOTO 178 c endif c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c...Generate synthetic spectrum from loaded interferograms (with U) c if(KK.eq.'U'.or.KK.eq.'u')then c if(bw.eq.0.d0)bw=1.d0 if(df.eq.0.d0)df=0.001d0 fs_st=freq(1)-0.5d0*BW fs_en=freq(nscans)+0.5d0*BW if(powamp.ne.'P'.and.powamp.ne.'A')powamp='P' if(nsfft.ge.0)nsfft=-4 c c...menu c dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) dummy=displaycursor($gcursoron) 59 CALL clearscreen($GCLEARSCREEN) c nsypt=1.d0+(fs_en-fs_st)/df WRITE(*,233) 233 FORMAT(1x/'_____S y n t h e t i c S p e c t r u m',40(1H_)//) WRITE(*,57)df,nsypt,natt,iabs(nsfft),powamp,bw 57 FORMAT(1x/ * ' 1: Synthetic point spacing (MHz) = ',F11.4,T60, * '[',i7' pts]'/ * ' 2: Decay attenuation (-ve) = ',I6/ * ' 3: n for FFT = ',I6/ * ' 4: Power/Amplitude = ',5x,a/ * ' 5: Coaddition width (MHz) = ',f11.4) if(intcut.eq.0)then write(*,234)intcut 234 format( * ' 6: Interf. cutoffs (auto/man =0/1) = ',I6) else write(*,235)intcut,nskipa,nskipb 235 format( * ' 6: Interf. cutoffs (auto/man =0/1) = ',I6,T60, * '[',2i7,']') endif write(*,236)iptick 236 format( * ' 7: Ticks at pump frequencies = ',I6// * ' 8: SYNTHESISE'//8x, * ' Select 1-8 (ENTER=Exit) ..... ',$) C C READ(*,'(I2)',err=59)I if(i.lt.0.or.i.gt.8)goto 59 c if(i.eq.1)then write(*,'(1x// * '' New point spacing for spectrum (MHz) = '',$)') read(*,'(F20.10)',err=59)freqv if(freqv.le.0.d0)goto 59 n=1.d0+(fs_en-fs_st)/freqv if(n.lt.10.or.n.gt.maxsyn)goto 59 df=freqv goto 59 endif c if(i.eq.2)then write(*,511) 511 format(1x//5x, * 'Interferograms can be multiplied by a decaying exponential ', * 'function to'/5x,'limit the effects of noise.'//5x, * 'The (-ve) decay constant for the exponential is to be ', * 'expressed as'/5x, * 'the number of interferogram points and the ', * 'default setting is a third'/5x, * 'of average interferogram length.'//5x, * 'New time-domain attenuation parameter = ',$) read(*,'(i10)',err=59)n if(n.ne.0.and.n.ge.-10)goto 59 natt=n goto 59 endif c c...n factor for FFT c if(i.eq.3)then write(*,'(1x// * '' New n for FFT (1...5) = '',$)') read(*,'(i10)',err=59)n if(n.lt.1.or.n.gt.5)goto 59 nsfft=-n goto 59 endif c c...power or amplitude c if(i.eq.4)then if(powamp.eq.'P')then powamp='A' else powamp='P' endif goto 59 endif c c...coaddition width c if(i.eq.5)then write(*,'(1x// * '' New Coaddition (full) bandwidth (MHz) = '',$)') read(*,'(F20.10)',err=59)freqv if(freqv.le.0.01d0.or.freqv.gt.1000.d0)goto 59 bw=freqv fs_st=freq(1)-0.5d0*BW fs_en=freq(nscans)+0.5d0*BW goto 59 endif c c...set interferogram limits c if(i.eq.6)then intcut=intcut+1 if(intcut.eq.2)then intcut=0 goto 59 endif 530 write(*,'(1x// * '' Interferogram cutoffs for synthesis = '',$)') read(*,*,err=530)nskipa,nskipb if(nskipa.le.0)goto 530 if(nskipb.gt.nrep)goto 530 if(nskipb-nskipa.lt.nrep/20)goto 530 goto 59 endif c c...ticks c if(i.eq.7)then iptick=iptick+1 if(iptick.ge.2)iptick=0 goto 59 endif c c...synthesise spectrum c if(i.eq.8)then if(nsypt.gt.maxsyn)then write(*,237)nsypt,maxsyn 237 format(1x//' **** ERROR: You need',i10, * ' points in the synthetic spectrum'/ * ' program limit is',i10,' points'// * ' Synthesis will be possible if you increase ', * 'frequency step'//) pause goto 59 endif c dummy=displaycursor($gcursoroff) do 500 n=1,nsypt synspe(n)=0.d0 500 continue c write(*,'(1x/)') open(2,file='ifpars.out',status='unknown') write(2,600) 600 format(80(1H-)/'Trace of automatic IF subtraction on spectrum', * ' synthesis'/80(1h-)// * 'Interferogram:'/' Niter CHIsq ALAMDA ', * ' ampl phase shift nfirst nlast'/) c do 501 n=1,nscans c write(*,'('' Interferogram'',i5,'': '',$)')n fcent=freq(n) jj=ipoint(n) nrep=nreps(jj) if(intcut.eq.0)then nskips=nskip(jj) nskipe=nskip1(jj) else nskips=nskipa nskipe=nrep-nskipb endif tstep=tsteps(jj) vstep=vsteps(jj) fifkhz=fikhzs(jj) fifmhz=0.001d0*fifkhz indxi=intpt(jj,1)-1 do 503 nn=1,nrep indxi=indxi+1 idata(nn)=intall(indxi) 503 continue c c___fit the IF parameters to use for subtraction c nauts=0.73*nrep nautf=0.97*nrep c sumint=0.d0 slow= 1.d+20 shigh=-1.d+20 do 1554 j=nauts,nautf sumint=sumint+idata(j) if(idata(j).lt.slow)slow=idata(j) if(idata(j).gt.shigh)shigh=idata(j) 1554 continue autoam=0.5d0*(shigh-slow) autoph=45.d0 autosh=sumint/(nautf-nauts+1) c macon=5 acon(1)=autoam acon(2)=autoph acon(3)=autosh acon(4)=0.d0 acon(5)=0.d0 lista(1)=1 lista(2)=2 lista(3)=3 mfit=3 nca=maxpar jjj=0 do 1557 j=nauts,nautf jjj=jjj+1 xpfit(jjj)=real(j)*tstep*1.d6 ypfit(jjj)=idata(j) sigpts(jjj)=1.d0 1557 continue npfit=nautf-nauts+1 c c c-----------Fitting loop-------\ c ntotit=0 1564 niter=0 alamda=-1.d0 CALL MRQMIN(XPFIT,YPFIT,SIGpts,NPFIT,Acon,MAcon,LISTA,MFIT, <----- * COVAR,ALPHA,NCA,CHISQ,SUBIF,ALAMDA) c 1556 OCHISQ=CHISQ niter=niter+1 ntotit=ntotit+1 CALL MRQMIN(XPFIT,YPFIT,SIGpts,NPFIT,Acon,MAcon,LISTA,MFIT, <----- * COVAR,ALPHA,NCA,CHISQ,SUBIF,ALAMDA) c c...renormalize amplitude and phase if necessary if(acon(1).lt.0.d0)then acon(1)=-acon(1) acon(2)=acon(2)+180.d0 endif if(acon(2).lt.0.d0)then acon(2)=acon(2)+(int(abs(acon(2)/360.d0))+1)*360.d0 endif if(acon(2).gt.360.d0)then acon(2)=acon(2)-(int(abs(acon(2)/360.d0)))*360.d0 endif c IF (ABS(OCHISQ-CHISQ).GT.1.D-1)GOTO 1556 c if(niter.gt.40)goto 1564 if(ALAMDA.gt.1.d0)then acon(1)=autoam goto 1564 endif c c-----------Fitting loop-------/ c rifamp=acon(1) rifph=acon(2) rifsh=acon(3) c WRITE(2,1558)n,ntotit,CHISQ,ALAMDA,rifamp,rifph,rifsh, * nauts,nautf 1558 FORMAT(I4,':',i4,1pE16.7,E9.1,0PF11.2,f9.2,f10.2,i9,i7) c c___subtract IF c riffr=2.d0*pi*fifmhz riffr=riffr*tstep*1.d6 phase=rifph*2.d0*pi/360.d0 do 1555 j=1,nrep rif=rifamp*dcos(riffr*real(j)+phase)+rifsh idata(j)=idata(j)-rif 1555 continue c c___correct decay time c if(natt.lt.0)then do 506 j=1,nrep freqv=DBLE(J)/DBLE(Natt) if(dabs(freqv).gt.10)freqv=dsign(10.d0,freqv) idata(j)=iDATA(j)*DEXP(freqv) 506 continue endif c c...do the FFT c nfft=nsfft CALL FFTEXE <----- fstepm=dsign(1.d0,fifmhz)*fstep*0.001d0 c if(powamp.eq.'A')then do 507 j=1,npts if(p(j).ne.0.d0)p(j)=sqrt(p(j)) 507 continue endif c c Number of synthetic point for frequency f: ns=1+(f-fs_st)/df c Frequency of synthetic point ns: f=fs_st+(ns-1)*df c Frequency of point in P: f=(fcent-FIFmhz)+(n-1)*fstepm c Number of point in P for frequency f: np=1+(f-fcent+FIFmhz)/fstepm c c NST,NEN are frequency limits in the synthetic spectrum to be updated c from the current spectrum in P. Intensity values from P are c derived by interpolation c c nst=(fcent-bw*0.5d0-fs_st)/df nen=(fcent+bw*0.5d0-fs_st)/df do 502 nn=nst,nen freqv=fs_st+(nn-1)*df ip=NINT(1.d0+(freqv-fcent+FIFmhz)/fstepm) fp=(fcent-FIFmhz)+(ip-1)*fstepm c c Linear interp. c if(ip.lt.npts.and.ip.gt.1)then c synspe(nn)=synspe(nn)+p(ip)+ c * (freqv-fp)*(p(ip+1)-p(ip))/fstepm c endif c c Quadratic interp. if(ip.lt.npts.and.ip.gt.1)then x0=fp-fstepm x1=fp x2=fp+fstepm CL0=(freqv-x1)*(freqv-x2)/((x0-x1)*(x0-x2)) CL1=(freqv-x0)*(freqv-x2)/((x1-x0)*(x1-x2)) CL2=(freqv-x0)*(freqv-x1)/((x2-x0)*(x2-x1)) synspe(nn)=synspe(nn)+CL0*P(ip-1)+CL1*P(ip)+CL2*P(ip+1) endif 502 continue c 501 continue write(2,'(1x/80(1H-))') close(2) c c...add ticks at pump frequencies (if specified) c if(iptick.eq.1)then synmax=-1.e+20 synmin=1.e+20 do 550 n=1,nsypt if(synspe(n).gt.synmax)synmax=synspe(n) if(synspe(n).lt.synmin)synmin=synspe(n) 550 continue do 551 n=1,nscans nnn=1+nint((freq(n)-fs_st)/df) if(nnn.gt.0.and.nnn.le.nsypt)then synspe(nnn)=synspe(nnn)-0.05*(synmax-synmin) endif 551 continue endif c c...output c open(8,file='synspe.out',status='unknown') write(8,509)df,bw,nskips,nrep-nskipe,natt,iabs(nsfft),powamp 509 format('SPECTRUM SYNTHESIS: df=',F5.3,' BW=',f4.2, * ' use[',i5,',',i5,'] decay=',i5,' n=',i1,','A1) freqv=fs_st-df do 508 n=1,nsypt freqv=freqv+df write(8,512)freqv,synspe(n) 512 format(f10.4,1PE12.4) 508 continue close(8) c write(*,'(1x//5x,''Synthetic spectrum has been written to:'' * /14x,'' '',$)') dummy=settextcolor(int2(12)) write(*,'(a)')dirnam(1:len_trim(dirnam))//'\synspe.out' dummy=settextcolor(int2(ntextc)) WRITE(*,'(5x,''which can be further processed with SVIEW_L.''/ * 5x,''A trace of IF subtraction has also been '', * ''written to file IFPARS.OUT.''/// * 5x,''Press ENTER to continue '',$)') 200 ik=inkey(ik) if(ik.ne.13)goto 200 c dummy=displaycursor($gcursoron) goto 59 endif c C...return to main display c dummy=displaycursor($gcursoroff) CALL clearscreen($GCLEARSCREEN) CALL settextposition(8,int2(itxt+1),curpos) call outtext(emplin) GOTO 178 c endif c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c IF(IK.NE.13)GOTO 77 C C...exit control C DUMMY= SETTEXTCOLOR(int2(ntextc)) dummy4 = setbkcolor(ntextb) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(emplin) DUMMY=SETTEXTCOLOR(int2(15)) dummy4=setbkcolor (12) WRITE(outstr,'(A)')' ARE YOU SURE (Y/N) ?' CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(outstr(1:21)) 916 IK=INKEY(N) KK=CHAR(IK) IF(KK.EQ.'Y'.OR.KK.EQ.'y')GOTO 915 IF(KK.NE.'N'.AND.KK.NE.'n')GOTO 916 C DUMMY=SETTEXTCOLOR(int2(ntextc)) dummy4 = setbkcolor(ntextb) CALL settextposition(1,int2(itxt+1),curpos) call outtext(lwork1) GOTO 77 c 915 dummy=setexitqq(qwin$exitnopersist) stop c RETURN END c c C_____________________________________________________________________________ C subroutine marsca(f1,f2,FCENT,FIF,tstep) c c Routine to plot and label the marker scale for frequency limits F1 - move to beginning/end of spectra c D,F - alternative keys to / c W,Z - change vertical scaling c 1-9 - goto displayed spectrum 1-9 c F1-F9 - display parameters of spectrum 1-9 c - reread input data c - go back to previous screen c 77 IK=INKEY(N) KK=CHAR(IK) c c...terminate program (with Q) c if(KK.eq.'Q'.or.kk.eq.'q')then DUMMY= SETTEXTCOLOR(int2(ntextc)) dummy4 = setbkcolor(ntextb) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(emplin) DUMMY=SETTEXTCOLOR(int2(15)) dummy4=setbkcolor (12) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(' ARE YOU SURE YOU WANT TO EXIT (Y/N) ?') 916 IK=INKEY(N) KK=CHAR(IK) IF(KK.EQ.'Y'.OR.KK.EQ.'y')GOTO 915 IF(KK.NE.'N'.AND.KK.NE.'n'.and.kk.ne.'Q'.and.kk.ne.'q')GOTO 916 C DUMMY=SETTEXTCOLOR(int2(ntextc)) dummy4 = setbkcolor(ntextb) CALL settextposition(1,int2(itxt+1),curpos) call outtext(lwork1) GOTO 179 c 915 dum2=setexitqq(qwin$exitnopersist) stop endif c c...go to the left in spectra (with A) c if(KK.eq.'A'.or.kk.eq.'a')then nmark=nmark-1 if(nmark.gt.nscans-5)nmark=nscans-5 if(KK.eq.'A')nmark=nmark-4 if(nmark.lt.5)nmark=5 goto 179 endif c c...go to the right in spectra (with S) c if(KK.eq.'S'.or.kk.eq.'s')then nmark=nmark+1 if(nmark.lt.6)nmark=6 if(KK.eq.'S')nmark=nmark+4 if(nmark.gt.nscans-4)nmark=nscans-4 goto 179 endif c c...first,last spectrum or D,F c if(ik.eq.-71.or.kk.eq.'D'.or.kk.eq.'d')then nmark=5 goto 179 endif c if(ik.eq.-79.or.kk.eq.'F'.or.kk.eq.'f')then nmark=nscans-4 goto 179 endif c c...go to FFT of selected interferogram c if(KK.eq.'1'.or.KK.eq.'!')then ncall=nstart goto 800 endif if(KK.eq.'2'.or.KK.eq.'@')then ncall=nstart+1 goto 800 endif if(KK.eq.'3'.or.KK.eq.'#')then ncall=nstart+2 goto 800 endif if(KK.eq.'4'.or.KK.eq.'$')then ncall=nstart+3 goto 800 endif if(KK.eq.'5'.or.KK.eq.'%')then ncall=nstart+4 goto 800 endif if(KK.eq.'6'.or.KK.eq.'^')then ncall=nstart+5 goto 800 endif if(KK.eq.'7'.or.KK.eq.'&')then ncall=nstart+6 goto 800 endif if(KK.eq.'8'.or.KK.eq.'*')then ncall=nstart+7 goto 800 endif if(KK.eq.'9'.or.KK.eq.'(')then ncall=nstart+8 goto 800 endif c goto 801 800 if(ncall.gt.nscans)goto 801 if(iarch.eq.0)then ncall1=ncall else ncall1=ipoint(ncall) endif nrep=nreps(ncall1) nskips=nskip(ncall1) nskipe=nskip1(ncall1) ncolch=nrep-nskipe vstep=vsteps(ncall1) tstep=tsteps(ncall1) iseen(ncall)=11 ndet=ndets(ncall1) fifkhz=fikhzs(ncall1) fifmhz=0.001d0*fifkhz call looksp(ncall,ncall1) <----- GOTO 179 801 continue c c...change of vertical scaling - zoom out (with Z) c if(KK.eq.'Z'.or.kk.eq.'z')then wmult=wmult*1.1 if(kk.eq.'Z')wmult=wmult*1.364 goto 178 endif c c...change of vertical scaling - zoom in (with W) c if(KK.eq.'W'.or.kk.eq.'w')then wmult=wmult/1.1 if(kk.eq.'W')wmult=wmult/1.364 goto 178 endif c c...Display recording parameters, (with Fn) function key c if(ik.ge.-67.and.ik.le.-59)then ncall=iabs(ik)-59 ncall=nstart+ncall filnam=fnams(ncall) if(iarch.eq.0)then ncall1=ncall else ncall1=ipoint(ncall) endif dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) CALL clearscreen($GCLEARSCREEN) write(*,'(1x)') c write(*,'(1x//17x,''Spectral file: '',$)') DUMmy2=settextcolor(int2(12)) write(*,'(a)')filnam DUMmy2=settextcolor(int2(ntextc)) write(*,'(1x//1x,a/)')comnt(ncall1) WRITE(*,10)' Sample: ',SAMPL(ncall1) WRITE(*,10)' Time/Date: ',TIMD(ncall1) write(*,11)' No of points: ',nreps(ncall1) write(*,11)' Points skipped at beginning: ',nskip(ncall1) write(*,11)' Points skipped at end: ',nskip1(ncall1) write(*,12)' Microwave frequency (MHz): ',freq(ncall) write(*,12)' X-spacing (microseconds): ',tsteps(ncall1)/1.E-6 write(*,12)' Y-spacing (Volts) : ',vsteps(ncall1) write(*,11)' Number of averages: ',nave(ncall1) write(*,'(1x)') write(*,12)' Intermediate frequency (kHz): ',fikhzs(ncall1) write(*,11)' Detection type: ',ndets(ncall1) c write(*,'(1x//30x,''Press E N T E R to continue'')') c 10 format(1x,2a) 11 format(1x,a,i7) 12 format(1x,a,f20.12) c 108 IK=INKEY(J) IF(IK.NE.13)GOTO 108 dummy4=setbkcolor(1) DUMmy2=settextcolor(int2(7)) GOTO 178 endif c c...Help screen (with H) c if(KK.eq.'H'.or.KK.eq.'h')then dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) CALL clearscreen($GCLEARSCREEN) WRITE(*,232) 232 FORMAT(1x/'_____', * 'I N T E R F E R O G R A M S screen commands', * 31(1H_)/// * 9x,' A/S - spectra left/right'/ * 9x,' - first/last spectrum'/ * 9x,' D/F - alternative keys to /'/ * 9x,' W/Z - Y-scale zoom-in/zoom-out'// * 9x,' caps on/off - fast/slow change in the above'/// * 9x,' 1 to 9 - go to spectrum 1 to 9'/ * 9x,' to - display parameters of spectrum 1 to 9'/ * 9x,' - reread input data'/ * 9x,' ENTER - exit to the SUMMARY screen'/ * 9x,' Q - quit the program'/// * 9x,' Press ENTER to exit this HELP screen'/) 107 IK=INKEY(J) IF(IK.NE.13)GOTO 107 dummy4=setbkcolor(1) DUMmy2=settextcolor(int2(7)) GOTO 178 endif c c...Go back to the M/P to reread input data (with ESC) c IF(IK.eq.27)then irefr=1 dummy4=setbkcolor(ntextb) dummy=settextcolor(int2(ntextc)) call clearscreen($GCLEARSCREEN) return ENDIF c IF(IK.NE.13)GOTO 77 C C...exit back to SUMMAR C RETURN END C C_____________________________________________________________________________ C subroutine startg(iconf) c C This routine uses QWIN graphics and techniques from the CLEANWIN programming C example for CVF6 to avoid the full-screen startup of standard graphics, C while preserving a simple frame. C Note the use of the WIN32 routines MoveWindow, UpdateWindow, GetWindowLong, C SetWindowLong, GetHWndQQ - their operation and parameter values are not C really understood! c USE DFLIB USE DFWIN c RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy INTEGER*2 maxx,maxy,linofs,mymode,myrows,mycols,ifc,ifm,ifhb,ifb, * idelta,ivista integer*4 dummy4,i character fntnam*30,line*80 COMMON /limits/wxy,maxx,maxy,linofs,curpos,ixy, * mymode,myrows,mycols COMMON /gsets/wc,win,ifc,ifm,ifhb,ifb,idelta type (windowconfig)wc type (qwinfo)win logical status C C C...Determine the type of Windows (using a WIN32 function). Official Windows C version numbers are (os.dwMajorVersion,os.dwMinorversion,...): C C Server 2008 6.0.6001 ME 4.90 C Vista 6.0.6000 NT 4 4.0.1381 C Server 2003 5.2 98 4.10 C XP 5.1.2600 95 4.00 C 2000 5.0.2195 c type (t_OSVERSIONINFO)os os.dwOSVersionInfoSize=sizeof(os) status=GetVersionEx( os ) c write(*,*)os.dwPlatformId,os.dwMajorVersion,os.dwMinorversion c ntsys=os.dwPlatformId-1 C C...An alternative way to determine the type of Windows, but this has the C the disadvantage of a flashing COMMAND PROMPT screen C c status=SYSTEMQQ('ver>opsys.ver') c open(3,file='opsys.ver',status='unknown') c read(3,'(A)')line c read(3,'(A)')line c if(line(1:7).eq.'Windows')then c ntsys=0 c else c ntsys=1 c if(line(1:20).eq.'Microsoft Windows XP')ntsys=2 c endif c close(3) c status=SYSTEMQQ('del opsys.ver') C c...set the principal window parameters, as hardcoded below and c specified in the V6.CFG file c wc.numtextcols=80 wc.numtextrows=30 c open(3,file='c:\fft\v6.cfg',status='old',err=12) 7 read(3,'(a)')line if(line(1:1).eq.'!')goto 7 read(line,5)wc.numxpixels read(3,'(a)')line read(line,5)wc.numypixels read(3,'(a)')line fntnam=line(36:65) 5 format(35x,i4) c wc.fontsize=QWIN$EXTENDFONT wc.numcolors=-1 wc.extendfontname=trim(fntnam)//char(0) wc.extendfontsize=-1 c wc.extendfontattributes=0 8 read(3,'(a)')line if(line(1:1).eq.'!')goto 9 read(line,5)iattr if(iattr.lt.1.or.iattr.gt.15)then write(*,10)iattr 10 format(1x//' Extended font attribute from V6.CFG is',i5, * ', which is illegal (1-15 allowed)'// * ' **** TRY AGAIN! *****'//) pause stop endif if(iattr.eq. 1)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_NORMAL if(iattr.eq. 2)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_UNDERLINE if(iattr.eq. 3)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_BOLD if(iattr.eq. 4)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_ITALIC c if(iattr.eq. 5)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_FIXED_PITCH if(iattr.eq. 6)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_VARIABLE_PITCH c if(iattr.eq. 7)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_FF_ROMAN if(iattr.eq. 8)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_FF_SWISS if(iattr.eq. 9)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_FF_MODERN if(iattr.eq.10)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_FF_SCRIPT if(iattr.eq.11)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_FF_DECORATIVE c if(iattr.eq.12)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_ANSI_CHARSET if(iattr.eq.13)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_DEFAULT_CHARSET if(iattr.eq.14)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_SYMBOL_CHARSET if(iattr.eq.15)wc.extendfontattributes=wc.extendfontattributes * +QWIN$EXTENDFONT_OEM_CHARSET goto 8 c 9 close(3) iconf=1 goto 11 C C...open default sized window for a 800x600 screen if the configuration file C cannot be opened C 12 wc.numxpixels=800 wc.numypixels=540 wc.extendfontname='Courier New' iconf=0 C C...Kill menu C 11 DO I = 7,1,-1 STATUS= DELETEMENUQQ(I, 0) END DO C C...Kill status bar C i = clickqq( QWIN$STATUS ) C C...Kill scroll bars and unwanted features (note that the title seems possible C only on the (killed) daughter window and not on the framewindow C i = GetWindowLong( GetHWndQQ(QWIN$FRAMEWINDOW), GWL_STYLE ) i = ior( iand( i, not(WS_THICKFRAME) ), WS_BORDER ) i = iand( i, not(WS_MAXIMIZEBOX) ) k = SetWindowLong( GetHWndQQ(QWIN$FRAMEWINDOW), GWL_STYLE, i ) C i = GetWindowLong( GetHWndQQ(0), GWL_STYLE ) i = ior(iand( i, not(WS_CAPTION.or.WS_SYSMENU.or.WS_THICKFRAME)), & & WS_BORDER) k = SetWindowLong( GetHWndQQ(0), GWL_STYLE, i ) c c...Position window - scrren coordinates of the top left corner are specified c ifxed= GetSystemMetrics(sm_cxfixedframe) x thickness of frame around nonzizable window ifyed= GetSystemMetrics(sm_cyfixedframe) y thickness of frame around nonzizable window c win.x = -ifxed c win.y = -ifyed win.x = 0 win.y = 0 c c...Correct sizing parameters to take account of removal of the menu bar and of c the caption area for the child window - these vary with system font size and c are augmented with emprirically established IDELTA fudge parameter C C C W98,ME small fonts: ifc=19 ifm=19 idelta=5 -> sum=43 C W2000, small=100% fonts: ifc=19 ifm=19 idelta=5 -> sum=43 C W2000, large=125% fonts: ifc=24 ifm=24 idelta=4 -> sum=52 C W2000, custom=132% fonts: ifc=25 ifm=25 idelta=3 -> sum=53 C WindowsXP,small fonts: ifc=26 ifm=20 idelta=2 -> sum=48 C WindowsXP, 125% fonts: ifc=32 ifm=25 idelta=2 32,25,1,1,3,3 C Vista, 125% fonts: ifc=25 ifm=25 25,25,1,1,3,3 C C ifc = GetSystemMetrics(sm_cycaption) ht of caption area ifm = GetSystemMetrics(sm_cymenu) ht of single line menu bar ifhb = GetSystemMetrics(sm_cxborder) width of border ifb = GetSystemMetrics(sm_cyborder) height of border c write(*,*)ifc,ifm,ifb,ifhb,ifxed,ifyed c pause c if(os.dwMajorVersion.eq.6)then ivista=10 else ivista=0 endif c idelta=4 if(ifc.eq.19.and.ifc.eq.19)idelta=5 if(ifc.eq.25.and.ifc.eq.25)idelta=3 if(ifc.ge.26)idelta=2 c win.w = wc.numxpixels-(2*ifhb+ivista) win.h = wc.numypixels-(ifc+ifm+idelta+ivista) c win.type=qwin$set dummy4 = setwsizeqq(qwin$framewindow,win) status = getwsizeqq(QWIN$FRAMEWINDOW,QWIN$SIZECURR, win) c wc.numtextcols=80 wc.numtextrows=30 wc.title=' 'C c status=setwindowconfig(wc) if(.not.status)status=setwindowconfig(wc) C C...Magical Windows incantations to make style set above real (without C these commands the active window does not expand to the size of the C program framewindow) C i = MoveWindow( GetHWndQQ(0), -1, -1, 0, 0, int4(.TRUE.)) call clearscreen($GCLEARSCREEN) status = UpdateWindow(GETHANDLEFRAMEQQ()) c C pixel limits on x and y axes (0,maxx), (0,maxy) c maxx=wc.numxpixels-1 maxy=wc.numypixels-1 myrows=wc.numtextrows mycols=wc.numtextcols linofs=nint(real(maxy)/real(myrows))+1 c c write(*,*)maxx,maxy c write(*,*)myrows,mycols c pause c return end C c--------------------------------------------------------------------------- c subroutine looksp(nmark,nmark1) c c NOTE that the frequency axes for the FFT, for actual plotting, and c as displayed in numerical values are all different: c c NDET=1 FFT plot minus IF displayed values c c IF=+20.000 MHz 0 +25 -> 0 +25 -> -20.000 -> [-20.000 +5.000] + FCENT c c IF=-20.278 MHz 0 +25 -> -25 0 -> +20.278 -> [ -4.722 +20.278] + FCENT c | | c FSTART FEND c c Also remember: c 1/ for negative IF the FFT axis is not shifted but reversed c 2/ the 'plot' scale is in kHz, as are FSTART and FEND c USE DFLIB use globdat C RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy c c fcent = excitation frequency (MHz) c fstep = point spacing in frequency (kHz) c p = spectral points from the FFT c npts = number of spectral points out of the FFT c fmax = frequency of the last spectral point (kHz) c fmark = current marker frequency (kHz) c fincr = frequency increment per horizontal pixel (kHz) c nmark = position of spectrum defined by FMARK in order of frequency c nmark1 = the original position (in order of input) of spectrum defined by c FMARK c logical*2 true real*8 pi,phase,rifph,rifamp,rifsh,riffr, * autoam,autoph,autosh PARAMETER (Nmaxpt=131072,maxsmo=1001,maxpar=10, * nivols=7,nsect=20,true=.true.) PARAMETER (ntextc=0, ntextb=7, nbordc=15, ncursc=14) PARAMETER (pi=3.141592653589793d0) c integer*4 dummy4,idata(maxpts),ioldat(maxpts),itemp(maxpts) real*8 fcent,fstep,fstart,fend,top,bottom,f,fmark,fincr, * flast,fmean,fmax,fchang,frange,FMIN,FPLUS, * fnewp,fnewo,fnewm,flastp,flasto,flastm,smin,smax, * cursle,cursri,rinter,ssum,fprint,flims,flime, * x(nsect),y(nsect),xpeak,rspec real*8 fifmhz,fifkhz,fifmho,fifkho,fikhzs(maxspe) real p(nmaxpt),spol(maxsmo) INTEGER*2 dummy,maxx,maxy,LINOFS,inkey,ipoint(maxspe), * mymode,myrows,mycols character kk,outstr*41,filarc*30 character*80 emplin,lwork1,lwork2,lwork3 real*8 freq(maxspe),YYSTEP,RRSMAL,YSHIFT,F1,F2,t1,t2, * x1,y1,x2,y2 character fnams(maxspe)*12 real*4 detvol(maxspe,2),volint(maxspe,nivols),nave(maxspe) integer*2 iseen(maxspe),ndets(maxspe) integer*2 nreps(maxspe),nskip(maxspe),nskip1(maxspe) real*4 tsteps(maxspe),vsteps(maxspe) c real*8 xpfit(maxpts),ypfit(maxpts),sigpts(maxpts),acon(maxpar), * COVAR(maxpar,maxpar),ALPHA(maxpar,maxpar),CHISQ, * ALAMDA,ochisq integer*4 lista(maxpar) real*8 SL_ampl,SL_tau,SL_freq,SL_phase,SL_shift character comnt(maxspe)*50,sampl(maxspe)*20,timd(maxspe)*20 c EXTERNAL SUBIF,SUBLIN common /scans/freq,wmult,ipoint,filarc COMMON /limits/wxy,maxx,maxy,LINOFS,curpos,ixy, * mymode,myrows,mycols common /specf/p,fstep,npts,NFFT,NCALL common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /scans2/fikhzs,detvol,volint,iseen,nscans,ndets common /scan3/nreps,nskip,nskip1,tsteps,vsteps,nave common /sfiles/fnams common /peak/x,y common /smooth/ioldat,itemp,spol COMMON /plotda/bottom,top,RRSMAL,YYSTEP,YSHIFT common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto COMMON /SCTEXT/COMNT,SAMPL,TIMD c indxi=intpt(nmark1,1)-1 do 100 i=1,nrep indxi=indxi+1 idata(i)=intall(indxi) 100 continue cursm=nrep/500+1 curbi=10.*cursm curba=9.*cursm c fcent=freq(nmark) nauts=0 nautf=0 rifamp=0.d0 new=1 c c...preserve a copy of the interferogram in IOLDAT c do 1099 j=1,nrep ioldat(j)=idata(j) 1099 continue c c...apply the default decay correction c nhalf=0 c nhalf=-real(nrep)/3.d0 c do 1507 j=1,nrep c f=DBLE(J)/DBLE(NHALF) c if(dabs(f).gt.10)f=dsign(10.d0,f) c idata(j)=iDATA(j)*DEXP(f) c1507 continue c c c...subtract the IF by visiting the code in the interferogram viewing c screen below, which is executed with i c ncall=0 if(ndet.eq.1)goto 1560 1507 continue c c...do the first FFT c dummy4 = setbkcolor( ntextb ) call clearscreen($GCLEARSCREEN) call FFTEXE <----- samfre=1./tstep fnyq=samfre/2. fmax=fnyq*1.E-3 c htcut=0.5 CURINC=1.0 WRITE(emplin,'(80(1H ))') c c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...Start up the graphics C C C...definition of graphics viewport in pixel coordinates C call setviewport(int2(2),int2(2*LINOFS-2), * int2(maxx-2),int2(maxy-2*LINOFS+2)) <-vvvv c c...Find Y-limits c 851 if(fifmhz.ne.0.d0)then fstart=fifkhz-500.d0 fend= fifkhz+500.d0 else fstart=0.d0 fend=fstart+500.d0 endif fmark=0.5d0*(fstart+fend) fstep=dsign(1.d0,fifmhz)*dabs(fstep) c c...FLIMS,FLIME are frequency limits for the frequency axis (in kHz and c without adding fcent) c FLIMS=0 for +ve and zero IF (window is 0,FMAX) c FLIMS=FMAX for -ve IF (window is -FMAX,0) c flims=0.5d0*(nint(dsign(1.d0,fifkhz))-1)*fmax flime= flims+fmax c smin=1.d+20 smax=-1.d+20 f=-fstep do 111 i=1,npts f=f+fstep if(fifmhz.ge.0.d0)then if(f.lt.fstart)goto 111 if(f.gt.fend)goto 112 else if(f.gt.fend)goto 111 if(f.lt.fstart)goto 112 endif if(p(i).lt.smin)smin=p(i) if(p(i).gt.smax)smax=p(i) 111 continue 112 smaxol=smax SMIN=0 c top= smax+0.10d0*(smax-smin) bottom=smin-0.10D0*(smax-smin) YYSTEP=(top-bottom)/500.d0 RRSMAL=bottom-13.d0*YYSTEP itxt=(mycols-80)/2 c c...definition of graphics window in floating point coordinates to be used c for plotting c c The graphics viewport is assigned to be two pixels narrower on each c side than the screen, and without the number of pixels corresponding c to two lines at the top and one line at the bottom c c...Plot c 699 dummy=setwindow(TRUE,fstart,RRSMAL,fend,top) <-wwww fincr=(fend-fstart)/maxx c dummy=setcolor(int2(nbordc)) dummy4 = setbkcolor( 1 ) call clearscreen($GVIEWPORT) DUMMY2=SETCOLOR(int2( 0 )) CALL moveto_w(fstart,bottom,wxy) dummy=lineto_w(fstart,top) dummy=lineto_w(fend,top) DUMMY2=SETCOLOR(int2( nbordc )) dummy=lineto_w(fend,bottom) dummy=lineto_w(fstart,bottom) nfirst=0 DO 6 I=1,npts f=(i-1)*fstep if(fifmhz.ge.0.d0)then if(f.lt.fstart)goto 6 if(f.gt.fend)goto 697 else if(f.gt.fend)goto 6 if(f.lt.fstart)goto 697 endif if(nfirst.eq.0)then nfirst=1 RSPEC=p(i) IF(RSPEC.gt.top)rspec=top if(rspec.lt.0.d0)rspec=0.d0 CALL moveto_w(f,dble(rspec),wxy) endif RSPEC=p(i) IF(RSPEC.gt.top)rspec=top if(rspec.lt.0.d0)rspec=0.d0 dummy=lineto_w(f,RSPEC) 6 CONTINUE C C...marker scale C 697 yshift=bottom if(ndet.eq.1)then f1=fcent+fstart*0.001d0-FIFMHZ f2=fcent+ fend *0.001D0-FIFMHZ call marsca(f1,f2,fcent,FIFMHZ,0.0d0) <----- endif if(ndet.eq.0)then f1=fstart*0.001d0 f2= fend*0.001d0 call marsca(f1,f2,0.d0,0.d0,0.0d0) <----- endif c c...f-f_pump=0 marker (if within display range) c if(FIFkhz.gt.fstart.and.FIFkhz.lt.FEND)then DUMMY2=SETCOLOR(int2( 11 )) F1=bottom+(top-bottom)*0.1d0 CALL moveto_w(FIFkhz,F1,wxy) dummy2=lineto_w(FIFkhz,bottom) dummy2=lineto_w(FIFkhz+10.d0*fincr,RRSMAL) dummy2=lineto_w(FIFkhz-10.d0*fincr,RRSMAL) dummy2=lineto_w(FIFkhz,bottom) endif c c...cursor c dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,top,wxy) dummy2=lineto_w(FMARK,bottom) dummy=setwritemode($GPSET) CALL setlinestyle(int2(#FFFF)) c c...information lines c 771 dummy4=setbkcolor(ntextb) dummy=settextcolor(int2(ntextc)) FMIN=FCENT-FMARK/1000.D0 FPLUS=FCENT+FMARK/1000.D0 YVAL=P( NINT(FMARK/FSTEP)+1 ) c if(ndet.eq.0)then write(lwork1,'(F7.1,''kHz --> f-'',F11.4,'', f+'',F11.4, * '' MHz'',20X,''Y:'',F9.2)')fmark,FMIN,FPLUS,yval endif c if(ndet.eq.1)then fprint=fcent-fifmhz+fmark*0.001d0 write(lwork1,'(F8.1,''kHz --> f='',F11.4,'' MHz'',15x, * 17X,''Y:'',F11.4)')fmark-fifkhz,fprint,yval endif c CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork1) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(EMPLIN) if(new.eq.1)then dummy=settextcolor(int2(2)) CALL settextposition(2,19,curpos) call outtext(comnt(ipoint(nmark))) dummy=settextcolor(int2(ntextc)) new=-1 endif c write(lwork3,822)nfft, smax-smin, fend-fstart 822 format(' n=',i1,F12.2,'<- Y_range X_range ->',f8.2,' kHz', * 22x,'H = help') CALL settextposition(myrows,int2(itxt+1),curpos) call outtext(lwork3(1:79)) DUMMY=SETTEXTCOLOR(int2(12)) CALL settextposition(myrows,int2(itxt+56),curpos) CALL outtext(fnams(nmark)) DUMMY=SETTEXTCOLOR(int2(ntextc)) c c...options loop c C K,L - scroll cursor c , - move cursor by quarter screen c A,S - scroll spectrum horizontally c Q,E - change horizontal scaling c W,Z - change vertical scaling C C I - display interferogram C CTRL I - display interferogram information C N - change the FFT zero-filling parameter c P - show the FFT points c B - subtract measured peak c R - return to initial settings c H - help screen c U - ASCII dump of current FFT to file F.DAT C (end is the higher of 0.5MHz and end of display window) C Y - ASCII dump of FFT in mV units (only after those have been C defined by a previous peek at the interferogram with I) C c O - determine frequency of peak nearest the cursor c 0 - change bisection range for peak measurement c 9 - take cursor frequency as line frequency c = - central frequency of Doppler doublet (from last two lines c measured with 'O') c CTRL O - measure nearest peak, then fit and subtract it in time domain c - exit back to the calling routine c 77 IK=INKEY(N) KK=CHAR(IK) c IF(KK.EQ.'K'.OR.KK.EQ.'k')GOTO 710 IF(KK.EQ.'L'.OR.KK.EQ.'l')GOTO 711 IF(KK.EQ.',')GOTO 750 IF(KK.EQ.'Q'.OR.KK.EQ.'q')GOTO 721 IF(KK.EQ.'E'.OR.KK.EQ.'e')GOTO 720 IF(KK.EQ.'R'.OR.KK.EQ.'r')GOTO 730 IF(KK.EQ.'A'.OR.KK.EQ.'a')GOTO 740 IF(KK.EQ.'S'.OR.KK.EQ.'s')GOTO 760 IF(KK.EQ.'Z'.OR.KK.EQ.'z')GOTO 810 IF(KK.EQ.'W'.OR.KK.EQ.'w')GOTO 820 IF(KK.EQ.'I'.OR.KK.EQ.'i')GOTO 830 IF(KK.EQ.'H'.OR.KK.EQ.'h')GOTO 840 IF(KK.EQ.'N'.OR.KK.EQ.'n')GOTO 850 IF(KK.EQ.'P'.OR.KK.EQ.'p')GOTO 1300 IF(KK.EQ.'B'.OR.KK.EQ.'b')GOTO 1900 IF(KK.EQ.'O'.OR.KK.EQ.'o')GOTO 1400 IF(KK.EQ.'0'.OR.KK.EQ.')')GOTO 1450 IF(KK.EQ.'9'.OR.KK.EQ.'(')GOTO 1440 IF(KK.EQ.'='.OR.KK.EQ.'-')GOTO 1430 IF(IK.EQ.15)goto 1460 with CTRL O if(ik.eq.-71)goto 7730 with home IF(IK.eq.9)GOTO 2100 with CTRL I C C...ASCII dump of current spectrum (with U) and in mV (with Y) C if(kk.eq.'U'.or.kk.eq.'u'.or.kk.eq.'Y'.or.kk.eq.'y')then open(7,file='f.dat',status='unknown') if(ndet.eq.0)fcut=max(500.,fend) if(ndet.eq.1)fcut=2.d0 write(7,'(1H!,50(1H-)/1H!/''! FFT of interferogram from: '',a/ * 1H!/1H!,50(1H-)/1H!)')fnams(nmark) write(7,'(''! Pump frequency ='',f15.6/ * ''! Maximum offset ='',f15.6/1H!/ * ''! Maximum intensity ='',f15.6/1H!)') * fcent,fcut/1000.,smax if(ndet.eq.0)then write(7,'(''! DeltaF power f- f+''/1H!)') endif if(ndet.eq.1)then write(7,'(''! Frequency Intensity''/1H!)') endif c pmult=1.d0 if(kk.eq.'Y'.or.kk.eq.'y')then pmax=0.d0 do 1800 i=1,npts if(p(i).gt.pmax)pmax=p(i) 1800 continue pmult=1000.d0*vstep*(maxi-mini)/pmax endif f=0.d0 c if(ndet.eq.0)then do 3 i=1,npts fmin=fcent-f/1000. fplus=fcent+f/1000. write(7,'(f8.5,1pe12.4,0P,2f12.4)')f/1000.,p(i)*pmult, * fmin,fplus f=f+fstep if(f.gt.fcut)goto 33 3 continue endif c if(ndet.eq.1)then if(fifmhz.gt.0.d0)then jj=1 jjj=npts jstep=1 endif if(fifmhz.lt.0.d0)then jj=npts jjj=1 jstep=-1 endif c DO 13 I=jj,jjj,jstep f=fcent+(i-1)*fstep*0.001d0-fifmhz if(f.lt.fcent-fcut)goto 13 if(f.gt.fcent+fcut)goto 13 write(7,'(f12.4,1pe12.4)')f,p(i)*pmult 13 CONTINUE endif c 33 write(7,'(1H!,50(1H-))') close(7) CALL settextposition(myrows/4,int2(itxt+1)+15,curpos) call outtext( * 'Spectrum from current FFT has been written to F.DAT') endif c IF(IK.NE.13)GOTO 77 C C...exit C 915 continue do 2201 j=1,nrep idata(j)=ioldat(j) 2201 continue return C C...Shift cursor to the left (with K) C 710 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,top,wxy) dummy2=lineto_w(FMARK,bottom) fmark=fmark-8.d0*fincr IF(KK.EQ.'k')fmark=fmark+7.d0*fincr IF(fmark.LT.fstart)fmark=fstart C 719 CALL moveto_w(FMARK,top,wxy) dummy2=lineto_w(FMARK,bottom) dummy=setwritemode($GPSET) CALL setlinestyle(int2(#FFFF)) GOTO 771 C C...Shift cursor to the right (with L) C 711 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,top,wxy) dummy2=lineto_w(FMARK,bottom) fmark=fmark+8.d0*fincr IF(KK.EQ.'l')fmark=fmark-7.d0*fincr IF(fmark.gT.fend)fmark=fend GOTO 719 C C...Center cursor, on second keypress move the cursor into the center of the C opposite screenhalf (with ,) C 750 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(FMARK,top,wxy) dummy2=lineto_w(FMARK,bottom) fmean=(fend+fstart)/2.d0 IF(fmark.EQ.Fmean)THEN IF(fLAST.LT.fmean)FMARK=FSTART+0.75*(fend-fstart) IF(fLAST.GE.fmean)FMARK=FSTART+0.25*(fend-fstart) fLAST=FMEAN GOTO 719 ENDIF fLAST=fMARK fmark=fmean GOTO 719 c c...zoom-in in frequency (with E) c 720 FRange=Fend-Fstart Fchang=0.10D0*FRange IF(KK.EQ.'e')Fchang=0.33d0*FRange Fstart=Fmark-Fchang Fend=Fmark+Fchang c 698 if(fstart.lt.flims)fstart=flims if(fend.gt.flime)fend=flime IF(FMARK.LT.fstart)FMARK=Fstart IF(FMARK.GT.fend)FMARK=Fend c goto 699 C C...zoom-out in frequency (with Q) C 721 FRange=Fend-Fstart Fchang=1.D0*FRange IF(KK.EQ.'q')Fchang=0.25d0*FRange Fstart=Fstart-Fchang Fend=Fend+fchang GOTO 698 c c...restore original settings (with R) c 730 if(fifmhz.ne.0.d0)then fstart=fifkhz-500.d0 fend= fifkhz+500.d0 else fstart=0.d0 fend=500.d0 endif fmark=0.5d0*(fstart+fend) smax=smaxol 802 top= smax+0.10d0*(smax-smin) bottom=SMIN-0.10D0*(smax-smin) YYSTEP=(top-bottom)/500.d0 RRSMAL=bottom-13.d0*YYSTEP goto 699 c c...go to low-frequency part of transform (with HOME) c 7730 if(fifmhz.ge.0.d0)then fstart=0.d0 fend=500.d0 else fend=flime fstart=flime-500.d0 endif fmark=0.5d0*(fend+fstart) goto 802 c c...shift window to the left (with A) c 740 FRange=fend-fstart Fchang=FRange*0.5D0 IF(KK.EQ.'a')Fchang=FRange*0.1D0 fstart=fstart-Fchang fend=fend-Fchang FMARK=FMARK-Fchang IF(fstart.LT.flims)THEN call BEEPQQ(2000,250) fstart=flims fend=fstart+fchang fmark=0.5d0*(fstart+fend) CALL settextposition(1,int2(itxt+1),curpos) 101 FORMAT(1X,A1,$) ENDIF GOTO 699 C C...shift of viewing window to the right (with S) C 760 FRange=fend-fstart Fchang=FRange*0.5D0 IF(KK.EQ.'s')Fchang=FRange*0.1D0 fstart=fstart+Fchang fend=fend+Fchang FMARK=FMARK+Fchang IF(fend.gt.flime)THEN call BEEPQQ(2000,250) fend=flime fstart=fend-fchang fmark=0.5d0*(fstart+fend) CALL settextposition(1,int2(itxt+1),curpos) ENDIF GOTO 699 C C...zoom-out in height (with Z) C 810 sMULT=2.D0 IF(KK.EQ.'z')sMULT=1.1D0 smax=SMIN+sMULT*(smax-SMIN) GOTO 802 C C...zoom-in in height (with W) C 820 sMULT=0.5D0 IF(KK.EQ.'w')sMULT=0.95D0 smax=SMIN+sMULT*(smax-SMIN) GOTO 802 c c...Show the FFT points (with P) C 1300 dummy=setcolor(int2(nbordc)) dely= (5.0/maxy)*smax delx=0.75*(5.0/maxx)*(fend-fstart) DO 1301 i=1,Npts if(p(i).gt.top)goto 1301 f=(i-1)*fstep if(fifmhz.ge.0.d0)then if(f.lt.fstart)goto 1301 if(f.gt.fend)goto 1302 else if(f.gt.fend)goto 1301 if(f.lt.fstart)goto 1302 endif X1=f-DELX Y1=p(i)+DELY X2=f+DELX Y2=p(i)-DELY dummy=ellipse_w($GFILLINTERIOR,X1,Y1,X2,Y2) 1301 CONTINUE 1302 GOTO 77 c c...Show the interferogram information (with CTRL I) C 2100 dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) CALL clearscreen($GCLEARSCREEN) write(*,'(1x)') c write(*,'(1x//17x,''Spectral file: '',$)') DUMmy2=settextcolor(int2(12)) write(*,'(a)')fnams(nmark) DUMmy2=settextcolor(int2(ntextc)) write(*,'(1x//1x,a/)')comnt(nmark1) WRITE(*,2110)' Sample: ',SAMPL(nmark1) WRITE(*,2110)' Time/Date: ',TIMD(nmark1) write(*,2111)' No of points: ',nreps(nmark1) write(*,2111)' Points skipped at beginning: ',nskip(nmark1) write(*,2111)' Points skipped at end: ',nskip1(nmark1) write(*,2112)' Microwave frequency (MHz): ',freq(nmark) write(*,2112) * ' X-spacing (microseconds): ',tsteps(nmark1)/1.E-6 write(*,2112)' Y-spacing (Volts) : ',vsteps(nmark1) write(*,2111)' Number of averages: ',nave(nmark1) write(*,'(1x)') write(*,2112)' Intermediate frequency (kHz): ',fikhzs(nmark1) write(*,2111)' Detection type: ',ndets(nmark1) c write(*,'(1x//30x,''Press E N T E R to continue'')') c 2110 format(1x,2a) 2111 format(1x,a,i7) 2112 format(1x,a,f20.12) c 2114 IK=INKEY(J) IF(IK.NE.13)GOTO 2114 dummy4=setbkcolor(1) DUMmy2=settextcolor(int2(7)) goto 699 c c...Subtract the currently measured peak (with B) c c This option attempts to subtract the most recent peak that was measured c with the O option by fitting its waveform in the time domain C 1900 if(ndet.eq.0)goto 77 goto 77 c c...Take current marker frequency as measurement of line frequency (with 9) c 1440 xpeak=fmark FLASTM=FNEWM FLASTP=FNEWP flasto=fnewo fNEWM=fcent-xpeak*0.001d0 FNEWP=FCENT+XPEAK*0.001d0 fnewo=xpeak c if(ndet.eq.0)then write(lwork3,2441)xpeak,'peak',FNEWM,FNEWP 2441 format(F7.1,'kHz <-',a,'-> ',f10.4,' - + ',f10.4) endif if(ndet.eq.1)then fprint=fcent-fifmhz+xpeak*0.001d0 write(lwork3,1441)xpeak-fifkhz,'peak',fprint 1441 format(F8.1,'kHz <-',a,'-> ',f10.4) endif c DUMMY=SETTEXTCOLOR(int2(1)) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork3) write(lwork3,'(''pump='',f10.4)')fcent DUMMY=SETTEXTCOLOR(int2(12)) CALL settextposition(1,int2(itxt+51),curpos) CALL outtext(lwork3(1:15)) c CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) DUMMY=SETTEXTCOLOR(int2(ntextc)) c DUMMY=SETCOLOR(int2(14)) y1=top y2=0.d0 call moveto_w(xpeak,DBLE(y1),wxy) dummy=lineto_w(xpeak,DBLE(y2)) dummy=setcolor(int2(15)) c goto 77 c c...Frequency of peak nearest the cursor (with O) C 1400 FLASTM=FNEWM FLASTP=FNEWP flasto=fnewo ypeakl=ypeak call pf(fmark,ypeak,xpeak,errorx,htcut,xbot,fhalm,fhalp) <----- fNEWM=fcent-xpeak*0.001d0 FNEWP=FCENT+XPEAK*0.001d0 fnewo=xpeak C C...the quantity Er.fit is error in frequency determined from straight line C fit in Hz - in practice if it exceeds 10 then appreciable curvature C is present in bisector points for some rectifying action to be taken C c if(ndet.eq.0)then write(lwork3,2410)xpeak,'peak',FNEWM,FNEWP,errorx*1000.,ypeak 2410 format(F7.1,'kHz <-',a,'-> ',f10.4,' - + ',f10.4, * ' (Er.fit=',f6.1,') Y:',f9.2) endif if(ndet.eq.1)then fprint=fcent-fifmhz+xpeak*0.001d0 write(lwork3,1410)xpeak-fifkhz,'peak',fprint,errorx*1000.,ypeak 1410 format(F8.1,'kHz <-',a,'-> ',f10.4, * ' (Er.fit=',f6.1,')',13x,' Y:',f11.4) endif c DUMMY=SETTEXTCOLOR(int2(1)) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork3) WRITE(lwork3,'(67x,''FWHH:'',F8.2)')abs(fhalp-fhalm) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(lwork3) DUMMY=SETTEXTCOLOR(int2(ntextc)) c c...bisection points dely= (5.0/maxy)*smax delx=0.75*(5.0/maxx)*(fend-fstart) DUMMY=SETCOLOR(int2(9)) DO 1401 i=1,Nsect X1=x(i)-DELX Y1=y(i)+DELY X2=x(i)+DELX Y2=y(i)-DELY if(y1.lt.top.and.y2.lt.top)then * dummy=ellipse_w($GFILLINTERIOR,X1,Y1,X2,Y2) 1401 CONTINUE c c...fitted line DUMMY=SETCOLOR(int2(12)) y2=0.4*ypeak if(y2.lt.top)then call moveto_w(xpeak,DBLE(ypeak),wxy) dummy=lineto_w(DBLE(xbot),DBLE(y2)) endif c c...vertical bar DUMMY=SETCOLOR(int2(14)) y1=ypeak y2=0.d0 if(y1.lt.top)then call moveto_w(xpeak,DBLE(y1),wxy) dummy=lineto_w(xpeak,DBLE(y2)) endif c c...horizontal bars at 0, 0.5 and 1.0 Ymax x1=xpeak-5.d0*delx x2=xpeak+5.d0*delx y1=ypeak if(y1.lt.top)then call moveto_w(DBLE(x1),DBLE(y1),wxy) dummy=lineto_w(DBLE(x2),DBLE(y1)) endif Y1=0. call moveto_w(DBLE(x1),DBLE(y1),wxy) dummy=lineto_w(DBLE(x2),DBLE(y1)) c c...FWHH bar y1=ypeak*0.5 if(y1.lt.top)then call moveto_w(DBLE(fhalm),DBLE(y1),wxy) dummy=lineto_w(DBLE(fhalp),DBLE(y1)) endif DUMMY=SETCOLOR(int2(15)) c goto 77 c c...Mean frequency of last two measured peaks C e.g. for Doppler doublets (with =) C Intensity weighted mean of last two measured peaks C e.g. for averaging over unassigned splittings (with -) C 1430 if(kk.eq.'=')then fmin= (flastm+fnewm)*0.5d0 fplus=(flastp+fnewp)*0.5d0 foffs=(flasto+fnewo)*0.5d0 else fmin= (flastm*sqrt(ypeakl)+fnewm*sqrt(ypeak))/ * (sqrt(ypeakl)+sqrt(ypeak)) fplus=(flastp*sqrt(ypeakl)+fnewp*sqrt(ypeak))/ * (sqrt(ypeakl)+sqrt(ypeak)) foffs=(flasto*sqrt(ypeakl)+fnewo*sqrt(ypeak))/ * (sqrt(ypeakl)+sqrt(ypeak)) endif c if(ndet.eq.0)then write(lwork3,2411)foffs,'mean',fmin,fplus,abs(flasto-fnewo) 2411 format(F7.1,'kHz <-',a,'-> ',f10.4,' - + ',f10.4,' splitting', * f7.2,' kHz') endif if(ndet.eq.1)then fprint=fcent-fifmhz+foffs*0.001d0 write(lwork3,1411)foffs-fifkhz,'mean',fprint,abs(flasto-fnewo) 1411 format(F8.1,'kHz <-',a,'-> ',f10.4,' splitting', * f7.2,' kHz') endif c DUMMY=SETTEXTCOLOR(int2(9)) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork3) DUMMY=SETTEXTCOLOR(int2(ntextc)) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) if(kk.eq.'-')then CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(' intensity weighted') endif goto 77 c c...Change bisection range for peak frequency determination (with 0) c 1450 CALL clearscreen($GCLEARSCREEN) 1451 write(*,1452)htcut 1452 format(3x/' Current height range for bisection:',f6.2/ * 25x,' New range: ',$) read(*,'(F7.4)',ERR=1451)ypeak if(ypeak.le.0.0.or.ypeak.ge.1.0)goto 1451 htcut=ypeak dummy=setwindow(TRUE,fstart,bottom, <-wwww * fend,top) call clearscreen($GCLEARSCREEN) GOTO 699 c c c...Measure peak frequency closest to the cursor, fit a corresponding c component in the time domain and subtract it from the interferogram c (with CTRL O) c 1460 call pf(fmark,ypeak,xpeak,errorx,htcut,xbot,fhalm,fhalp) <----- c fNEWM=fcent-dble(xpeak)*0.001d0 c FNEWP=FCENT+dble(XPEAK)*0.001d0 c maxypf=-1000000000 minypf= 1000000000 jjj=0 do 1461 j=nskips+1,nrep-nskipe jjj=jjj+1 xpfit(jjj)=real(jjj)*tstep*1.d6 ypfit(jjj)=idata(j) sigpts(jjj)=1.d0 if(idata(j).gt.maxypf)maxypf=idata(j) if(idata(j).lt.minypf)minypf=idata(j) 1461 continue npfit=nrep-nskipe-nskips taumax=npfit*tstep*1.d6 tstart=nskips*tstep*1.d6 autoam=0.5d0*(maxypf-minypf) c SL_ampl=autoam SL_tau=0.3333d0*taumax SL_freq=XPEAK*0.001d0 SL_phase=45.d0 SL_shift=0.d0 c macon=5 acon(1)=SL_ampl acon(2)=SL_tau acon(3)=SL_freq acon(4)=SL_phase acon(5)=SL_shift mfit=3 lista(1)=1 lista(2)=2 lista(3)=4 nca=maxpar c write(*,'(1x/)') write(*,*)' Y-axis limits =',minypf,maxypf write(*,*)'Measured line freq. XPEAK/kHz =',xpeak write(*,*)'Measured line freq. SL_freq/MHz =',SL_freq write(*,*)'Taumax, Tstart, Autoam: ',taumax,tstart,autoam c pause c write(*,1462) 1462 format(1x/' Iter. Chisq ALAMDA ', * 'Amplitude Tau Phase'/) c c Parameters in the call to the fitting routine: C C XPFIT(NPFIT), YPFIT(NPFIT) - data points to be fitted, parent vectors C can be longer than NPFIT C SIGpts(NPFIT) - standard deviations on YPFIT points C ACON(MACON) - parameters to be used in the fit, parent C vector can be longer than MACON c LISTA(MFIT) - pointers to parameters in ACON to be fitted, C parent vector can be longer than MFIT c COVAR(NCA,NCA) - work arrays of size NCA, only the MFITxMFIT c ALPHA(NCA,NCA) subblock is used, final call with ALAMDA=0 c returns covariance matrix in COVAR and c Hessian in ALPHA c CHISQ - chi squared for the fit c SUBIF - user supplied function for calculating Y for c a given X c ALAMDA - progress of fit parameter c c c------Line Fitting loop-------\ c ntotit=0 1464 niter=0 alamda=-1.d0 CALL MRQMIN(XPFIT,YPFIT,SIGpts,NPFIT,Acon,MAcon,LISTA,MFIT, <----- * COVAR,ALPHA,NCA,CHISQ,SUBLIN,ALAMDA) WRITE(*,1468)ntotit,CHISQ,ALAMDA,acon(1),acon(2),acon(4),acon(3) 1468 FORMAT(i3,1PE19.10,E11.2,0pf12.3,2f10.3,f12.5) c 1466 OCHISQ=CHISQ niter=niter+1 ntotit=ntotit+1 CALL MRQMIN(XPFIT,YPFIT,SIGpts,NPFIT,Acon,MAcon,LISTA,MFIT, <----- * COVAR,ALPHA,NCA,CHISQ,SUBLIN,ALAMDA) c c...renormalize amplitude and phase if necessary if(acon(1).lt.0.d0)then acon(1)=-acon(1) acon(4)=acon(4)+180.d0 endif if(acon(4).lt.0.d0)then acon(4)=acon(4)+(int(abs(acon(4)/360.d0))+1)*360.d0 endif if(acon(4).gt.360.d0)then acon(4)=acon(4)-(int(abs(acon(4)/360.d0)))*360.d0 endif WRITE(*,1468)ntotit,CHISQ,ALAMDA,acon(1),acon(2),acon(4),acon(3) if(acon(2).gt.taumax)acon(2)=0.5*taumax if(acon(2).lt.0.05*taumax)acon(2)=0.5*taumax if(ntotit.eq.15)pause c IF (ABS(OCHISQ-CHISQ).GT.1.D+3)GOTO 1466 c c if(niter.gt.40)goto 1464 if(ALAMDA.gt.1.d0.and.ntotit.lt.100)then autoam=0.9*autoam acon(1)=autoam goto 1464 endif c c------Line Fitting loop-------/ c SL_ampl=acon(1) SL_tau=acon(2) SL_phase=acon(4) c write(*,'(1x/'' Press ENTER to continue '',$)') 201 ik=inkey(ik) if(ik.ne.13)goto 201 dummy=displaycursor($gcursoroff) jj=6 c c...subtract and exit c if(jj.eq.6)then riffr=2.d0*pi*SL_freq riffr=riffr*tstep*1.d6 phase=SL_phase*2.d0*pi/360.d0 do 1467 j=1,nrep rif=SL_ampl*dexp(-(j*tstep*1.d6-tstart)/SL_tau)* * dcos(riffr*real(j)+phase)+SL_shift idata(j)=idata(j)-rif 1467 continue c nfft=-iabs(nfft) call FFTEXE <----- goto 851 endif c goto 851 c c . . . . . . . . . . . . . . . . . . . . . . . . c c...display the interferogram (with I) and allow selection of discarded c points with scrollable cursors c 830 ICHANG=0 rmult=0.0 831 mini=1000000000 maxi=-1000000000 do 832 j=1,nrep if(j.lt.nskips)goto 832 if(j.gt.nrep-nskipe)goto 833 if(idata(j).lt.mini)mini=idata(j) if(idata(j).gt.maxi)maxi=idata(j) 832 continue c 833 if(maxi.le.mini)then dummy=setvideomode($DEFAULTMODE) write(*,'(1x/'' ---> ERROR: mini,maxi ='',2i12//)')mini,maxi stop endif toplim=1.05d0*(real(maxi)-real(mini))+dble(mini) botlim=dble(maxi)-1.05d0*(real(maxi)-real(mini)) dummy=setcolor(int2(nbordc)) dummy4 = setbkcolor( 1 ) dummy=setwindow(TRUE,1.0d0,dble(botlim),DBLE(nrep),dble(toplim)) <-wwww call clearscreen($GVIEWPORT) DUMMY2=SETCOLOR(int2( 0 )) CALL moveto_w(1.d0,dble(botlim),wxy) dummy=lineto_w(1.d0,dble(toplim)) dummy=lineto_w(dble(nrep),dble(toplim)) DUMMY2=SETCOLOR(int2( nbordc )) dummy=lineto_w(dble(nrep),dble(botlim)) dummy=lineto_w(1.d0,dble(botlim)) c c...part of interferogram skipped at the beginning c rinter=idata(1) if(idata(1).gt.maxi)rinter=maxi if(idata(1).lt.mini)rinter=mini CALL moveto_w(DBLE(1),rinter,wxy) c dummy=setcolor(int2(12)) if(nskips.gt.0)then DO 61 I=1,nskips RINTER=DBLE(Idata(I)) if(idata(i).gt.maxi)rinter=maxi if(idata(i).lt.mini)rinter=mini dummy=lineto_w(DBLE(I),RINTER) 61 CONTINUE endif c c...part of interferogram used for the FFT c dummy=setcolor(int2(15)) ssum=0.d0 DO 86 I=nskips+1,nrep-nskipe Rinter=DBLE(Idata(I)) dummy=lineto_w(DBLE(I),Rinter) ssum=ssum+rinter 86 CONTINUE c c...part of interferogram skipped at the end c dummy=setcolor(int2(12)) DO 62 I=nrep-nskipe+1,nrep Rinter=DBLE(Idata(I)) if(idata(i).gt.maxi)rinter=maxi if(idata(i).lt.mini)rinter=mini dummy=lineto_w(DBLE(I),Rinter) 62 CONTINUE issum=(nrep-nskipe-nskips-1) if(issum.eq.0)issum=1 ssum=ssum/issum dummy=setcolor(int2(9)) CALL moveto_w(dble(nskips+1),ssum,wxy) dummy=lineto_w(dble(nrep-nskipe),ssum) dummy=setcolor(int2(15)) c c...marker scale (first set values that are carried over in /plotda/ and c then set them back to those needed for frequency display) c YYSTEP=1.1d0*(maxi-mini)/500.d0 RRSMAL=botlim bottom=RRSMAL+13.d0*YYSTEP yshift=bottom c T1=1*tstep*1.D6 T2=nrep*tstep*1.D6 call marsca(t1,t2,0.0d0,0.0d0,dble(tstep)*1.D6) <----- c bottom=smin-0.10D0*(smax-smin) yshift=bottom YYSTEP=(top-bottom)/500.d0 RRSMAL=bottom-13.d0*YYSTEP c c...cursors c CURSLE=DBLE(NSKIPS) CURSRI=DBLE(NREP-NSKIPE) dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(cursle,dble(maxi),wxy) dummy2=lineto_w(cursle,dble(mini)) CALL moveto_w(cursri,dble(maxi),wxy) dummy2=lineto_w(cursri,dble(mini)) dummy=setwritemode($GPSET) CALL setlinestyle(int2(#FFFF)) C c...descriptive lines c dummy4=setbkcolor(ntextb) dummy=settextcolor(int2(ntextc)) WRITE(LWORK2,'(2A)') * ' A,S <-cursors-> K,L + - B N R U', * ' O H = Help' CALL settextposition(1,int2(itxt+1),curpos) call outtext(lwork2) CALL settextposition(2,int2(itxt+1),curpos) call outtext(emplin(1:79)) c write(lwork3,82) (maxi-mini)*vstep*1000, * nint(cursle),nint(cursri), * nint(cursle)*tstep*1.d6,nint(cursri)*tstep*1.d6, * (nint(cursri)-nint(cursle))*tstep*1.d6 82 format(' Y_used:',f9.2,'mV',19x, * 'X_used:',i5,'-',i5,' ->',F5.1,'-',F5.1,' =',F5.1,'us') c CALL settextposition(myrows,int2(itxt+1),curpos) call outtext(lwork3(1:79)) c c...options: - scrolling of left and right cursor K,L A,S C - change increment for scrolling the cursor + - C - background subtraction B C - compensation for rotational relaxation N C - return to original interferogram R c - ASCII dump of current interferogram c (to file T.DAT) U c - subtraction of background interferogr. O c - changing subtraction coefficient up/down arrow c 834 IK=INKEY(N) KK=CHAR(IK) IF(KK.EQ.'K'.OR.KK.EQ.'k')GOTO 870 IF(KK.EQ.'L'.OR.KK.EQ.'l')GOTO 875 IF(KK.EQ.'A'.OR.KK.EQ.'a')GOTO 880 IF(KK.EQ.'S'.OR.KK.EQ.'s')GOTO 885 IF(KK.EQ.'B'.OR.KK.EQ.'b')GOTO 1100 IF(KK.EQ.'R'.OR.KK.EQ.'r')GOTO 1200 IF(KK.EQ.'N'.OR.KK.EQ.'n')GOTO 1500 IF(KK.EQ.'I'.OR.KK.EQ.'i')GOTO 1560 IF(KK.EQ.'H'.OR.KK.EQ.'h')GOTO 1600 IF(KK.EQ.'O'.OR.KK.EQ.'o')GOTO 1700 if(ik.eq.-72.or.ik.eq.-80)goto 1750 IF(KK.EQ.'+'.OR.KK.EQ.'=')CURINC=CURINC*5.0 IF(KK.EQ.'-'.OR.KK.EQ.'_')CURINC=CURINC*0.2 c c...ASCII dump of current interferogram (with U) c if(kk.eq.'U'.or.kk.eq.'u')then open(7,file='t.dat',status='unknown') write(7,'(1H!,50(1H-)/1H!/''! Interferogram from file: '',a/ * 1H!/1H!,50(1H-)/1H!)')fnams(nmark) write(7,'(''! Pump frequency/MHz ='',f15.6)')fcent write(7,'(''! IF /MHz ='',f15.6)')fifmhz write(7,'(''! Detection type ='', i8 )')ndet write(7,'(''! Point spacing/microsec ='',f15.6)')tstep/1.e-6 write(7,'(''! Cutoff points/microsec ='',F15.6,'','',F15.6/ * ''! Last point/microsec ='',F15.6)') * nskips*tstep*1.e+6,(nrep-nskipe)*tstep*1.e+6, * nrep*tstep*1.e+6 write(7,'(''! Intensity limits/mV ='',F15.6,'','',F15.6)') * mini*1000.d0*vstep,maxi*1000.d0*vstep write(7,'(''!'')') write(7,'(''! X-limits (internal) ='',i8,7x,'','',i8)') * nskips,(nrep-nskipe) write(7,'(''! Y-limits (internal) ='',i8,7x,'','',i8)') * mini,maxi write(7,1802) 1802 format('!'/'! Time/us Y/mV Y_internal'/ * '! / \'/ * '! for FFT excluded'/ * '! | |') do 14 j=1,nrep if(j.lt.nskips.or.j.gt.nrep-nskipe) * write(7,115)real(j)*tstep/1.E-6, * real(idata(j))*1000.d0*vstep if(j.gt.nskips.and.j.lt.nrep-nskipe) * write(7,15)real(j)*tstep/1.E-6, * real(idata(j))*1000.d0*vstep,idata(j) if(j.eq.nskips.or.j.eq.nrep-nskipe) * write(7,116)real(j)*tstep/1.E-6, * real(idata(j))*1000.d0*vstep, * real(idata(j))*1000.d0*vstep,idata(j) 14 continue 115 format(f6.2,' * ',f10.3,11x,'*') 116 format(f6.2,2f10.3,i12) 15 format(f6.2,f10.3,' * ',i12) write(7,'(1H!,50(1H-))') close(7) CALL settextposition(myrows/4,int2(itxt+1)+15,curpos) call outtext( * 'Current interferogram has been written to T.DAT') endif IF(IK.NE.13)GOTO 834 c c...recalculate FFT if any changes in discarded points c 890 if(ichang.eq.1)THEN NFFT=-NFFT nskips=cursle nskipe=nrep-cursri CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(' --- R e c a l c u l a t i n g F F T --- ') CALL FFTEXE <----- goto 851 ELSE goto 699 ENDIF c C...Shift RIGHT cursor to the left (with K) C 870 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(cursri,dble(maxi),wxy) dummy2=lineto_w(cursri,dble(mini)) CURSRI=CURSRI-curbi*CURINC IF(KK.EQ.'k')CURSRI=CURSRI+curba*CURINC IF(CURSRI.LE.CURSLE)CURSRI=CURSLE+CURINC C 872 CALL moveto_w(cursri,dble(maxi),wxy) dummy2=lineto_w(cursri,dble(mini)) dummy=setwritemode($GPSET) CALL setlinestyle(int2(#FFFF)) ICHANG=1 CALL settextposition(myrows,int2(itxt+39),curpos) WRITE(OUTSTR,1550)nint(cursle),nint(cursri), * nint(cursle)*tstep*1.d6,nint(cursri)*tstep*1.d6, * (nint(cursri)-nint(cursle))*tstep*1.d6 1550 format('X_used:',i5,'-',i5,' ->',F5.1,'-',F5.1,' =',F5.1,'us') call outtext(outstr) GOTO 834 C C...Shift RIGHT cursor to the right (with L) C 875 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(cursri,dble(maxi),wxy) dummy2=lineto_w(cursri,dble(mini)) CURSRI=CURSRI+curbi*CURINC IF(KK.EQ.'l')CURSRI=CURSRI-curba*CURINC IF(CURSRI.GT.dble(NREP))CURSRI=dble(NREP) GOTO 872 c C...Shift LEFT cursor to the left (with A) C 880 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(cursle,dble(maxi),wxy) dummy2=lineto_w(cursle,dble(mini)) CURSLE=CURSLE-curbi*CURINC IF(KK.EQ.'a')CURSLE=CURSLE+curba*CURINC IF(CURSLE.LT.0.d0)CURSLE=0.d0 C 882 CALL moveto_w(cursle,dble(maxi),wxy) dummy2=lineto_w(cursle,dble(mini)) dummy=setwritemode($GPSET) CALL setlinestyle(int2(#FFFF)) ICHANG=1 CALL settextposition(myrows,int2(itxt+39),curpos) WRITE(OUTSTR,1550)nint(cursle),nint(cursri), * nint(cursle)*tstep*1.d6,nint(cursri)*tstep*1.d6, * (nint(cursri)-nint(cursle))*tstep*1.d6 call outtext(outstr) GOTO 834 c C...Shift LEFT cursor to the right (with S) C 885 dummy=setwritemode($GXOR) DUMMY2=SETCOLOR(int2( ncursc )) CALL setlinestyle(int2(#3333)) CALL moveto_w(cursle,dble(maxi),wxy) dummy2=lineto_w(cursle,dble(mini)) cursle=cursle+curbi*CURINC IF(KK.EQ.'s')cursle=cursle-curba*CURINC IF(cursle.ge.cursri)cursle=cursri-1.d0*CURINC GOTO 882 C C C...Subtraction of background by triple smoothing of interferogram with c least squares smoothing interval (with B) c 1100 CALL clearscreen($GCLEARSCREEN) 1101 WRITE(*,1102)nrep 1102 FORMAT(1x/'_____B A C K G R O U N D S U B T R A C T I O N', * 32(1H_)/// * 5x,'This option is equivalent to a high-pass filter and allows ' * 'removal'/ * 5x'of unwanted low-frequency fluctuations. The baseline is ', * 'defined '/ * 5x,'by triple least-squares smoothing.'// * 5x,'If the decay rate is also to be modified this should be ', * 'done only'/5x, * 'after background subtraction on the original interferogram.'/// * 5x,'Number of points in the interferogram ',18(1h.),i5/5x, * 'Number of points (>3 and odd) in smoothing interval .... ',$) READ(*,'(i5)',ERR=1100)NSPT IF(NSPT.LE.3.OR.NSPT.GT.maxsmo)GOTO 1100 IF((NSPT/2)*2.EQ.NSPT)GOTO 1100 WRITE(*,'(1X//'' S U B T R A C T I N G''//)') C do 1103 j=1,nrep ioldat(j)=idata(j) 1103 continue call baksub(nspt) <----- c ICHANG=1 GOTO 831 C C C...Compensate for rotational relaxation (with N) C 1500 CALL clearscreen($GCLEARSCREEN) 1501 WRITE(*,1502)nrep,nhalf 1502 FORMAT(1x/'_____MODIFY APPARENT RELAXATION TIME',44(1H_)/// * 5x,'This option multiplies the interferogram by an additional', * ' exponential'/ * 5x,'decay factor. Its decay half-time is to be specified as', * ' the number'/5x,'of interferogram points:'// * 5x,'- positive values amplify the interferogram tail, improving', * ' resolution,'/ * 5x,'- negative values attenuate the interferogram tail, ' * 'improving S/N.'/ * 5x,'- zero uses the original intrerferogram'/// * 12x,'Number of points in the interferogram ...',i6/ * 12x,'Number of points for halfdecay ..........',i6// * 12x,'New halfdecay points .................... ',$) READ(*,'(i6)',ERR=1500)N IF(N.eq.0)then nhalf=0 GOTO 1200 ENDIF IF(N.LE.10.AND.N.GE.-10)GOTO 1500 NHALF=N C do 1506 j=1,nrep f=DBLE(J)/DBLE(NHALF) if(dabs(f).gt.10)f=dsign(10.d0,f) idata(j)=iDATA(j)*DEXP(f) 1506 continue C ICHANG=1 GOTO 831 c c...Restore original interferogram (with R) c 1200 indxi=intpt(nmark1,1)-1 do 1201 j=1,nrep indxi=indxi+1 idata(j)=intall(indxi) 1201 continue ichang=1 goto 831 C C C...subtract IF: manual control (with I) C fast track (with i) C 1560 if(kk.eq.'i'.or.ncall.eq.0)goto 1561 dummy=displaycursor($gcursoron) c 1549 CALL clearscreen($GCLEARSCREEN) 1551 WRITE(*,1552) 1552 FORMAT(1x/'_____SUBTRACT THE INTERMEDIATE FREQUENCY',40(1H_)/// * 5x,'This option allows subtraction from the interferogram of ', * 'remnant'/ * 5x,'intermediate frequency signal at 20 MHz'//) c 1561 if(nauts.eq.0.and.nautf.eq.0)then nauts=0.73*nrep nautf=0.97*nrep endif c sumint=0.d0 slow= 1.d+10 shigh=-1.d+10 do 1554 j=nauts,nautf sumint=sumint+idata(j) if(idata(j).lt.slow)slow=idata(j) if(idata(j).gt.shigh)shigh=idata(j) 1554 continue autoam=0.5d0*(shigh-slow) autoph=45.d0 autosh=sumint/(nautf-nauts+1) if(jj.eq.1.or.rifamp.eq.0.d0.or.KK.eq.'i')then rifamp=autoam rifph=autoph rifsh=autosh endif if(KK.eq.'i'.or.ncall.eq.0)goto 1562 c WRITE(*,1553)nauts,nautf,rifamp,autoam,rifph,autoph,rifsh,autosh 1553 FORMAT( * ' 1: Points for AUTO pars. = ',I9,6x,i9/ * ' 2: Refit parameters (no subtraction)' // * 35x,'Current parameter AUTO hint'/ * ' 3: IF amplitude (in A*cos) = ',2F15.5/ * ' 4: IF phase (degrees) = ',2F15.5/ * ' 5: IF shift = ',2F15.5// * ' 6: Subtract and display'/// * ' Select 1-6 (ENTER=exit) ..... ',$) C READ(*,'(I2)',err=1549)jj if(jj.lt.0.or.jj.gt.6)goto 1549 c c...change auto cutoffs c if(jj.eq.1)then write(*,'(1x// * 6x,'' New auto cutoffs = '',$)') read(*,*,err=1549)j,jjj if(j.lt.1)goto 1549 if(jjj.gt.nrep)goto 1549 if(j.ge.jjj)goto 1549 nauts=j nautf=jjj goto 1549 endif c c...AUTO fit c 1562 if(jj.eq.2.or.KK.eq.'i'.or.ncall.eq.0)then c macon=5 acon(1)=rifamp acon(2)=rifph acon(3)=rifsh acon(4)=0.d0 acon(5)=0.d0 lista(1)=1 lista(2)=2 lista(3)=3 mfit=3 nca=maxpar jjj=0 do 1557 j=nauts,nautf jjj=jjj+1 xpfit(jjj)=real(j)*tstep*1.d6 ypfit(jjj)=idata(j) sigpts(jjj)=1.d0 1557 continue npfit=nautf-nauts+1 c if(KK.eq.'I')write(*,1559) 1559 format(1x/' Iter. Chisq ALAMDA ', * 'Amplitude Phase Yshift'/) c c Parameters in the call to the fitting routine: C C XPFIT(NPFIT), YPFIT(NPFIT) - data points to be fitted, parent vectors C can be longer than NPFIT C SIGpts(NPFIT) - standard deviations on YPFIT points C ACON(MACON) - parameters to be used in the fit, parent C vector can be longer than MACON c LISTA(MFIT) - pointers to parameters in ACON to be fitted, C parent vector can be longer than MFIT c COVAR(NCA,NCA) - work arrays of size NCA, only the MFITxMFIT c ALPHA(NCA,NCA) subblock is used, final call with ALAMDA=0 c returns covariance matrix in COVAR and c Hessian in ALPHA c CHISQ - chi squared for the fit c SUBIF - user supplied function for calculating Y for c a given X c ALAMDA - progress of fit parameter c c c--------IF Fitting loop-------\ c 1564 niter=0 alamda=-1.d0 CALL MRQMIN(XPFIT,YPFIT,SIGpts,NPFIT,Acon,MAcon,LISTA,MFIT, <----- * COVAR,ALPHA,NCA,CHISQ,SUBIF,ALAMDA) if(KK.eq.'I')WRITE(*,1558)niter,CHISQ,ALAMDA, * acon(1),acon(2),acon(3) 1558 FORMAT(i3,1PE19.10,E11.2,0p3f12.3) c 1556 OCHISQ=CHISQ niter=niter+1 CALL MRQMIN(XPFIT,YPFIT,SIGpts,NPFIT,Acon,MAcon,LISTA,MFIT, <----- * COVAR,ALPHA,NCA,CHISQ,SUBIF,ALAMDA) c c...renormalize amplitude and phase if necessary if(acon(1).lt.0.d0)then acon(1)=-acon(1) acon(2)=acon(2)+180.d0 endif if(acon(2).lt.0.d0)then acon(2)=acon(2)+(int(abs(acon(2)/360.d0))+1)*360.d0 endif if(acon(2).gt.360.d0)then acon(2)=acon(2)-(int(abs(acon(2)/360.d0)))*360.d0 endif if(KK.eq.'I')WRITE(*,1558)niter,CHISQ,ALAMDA, * acon(1),acon(2),acon(3) c IF (ABS(OCHISQ-CHISQ).GT.1.D-1)GOTO 1556 c if((KK.eq.'i'.or.ncall.eq.0).and.niter.gt.40)goto 1564 if(ALAMDA.gt.1.d0)then acon(1)=autoam goto 1564 endif c c--------IF Fitting loop-------/ c rifamp=acon(1) rifph=acon(2) rifsh=acon(3) autoam=acon(1) autoph=acon(2) autosh=acon(3) c if(KK.eq.'I')then write(*,'(1x/'' Press ENTER to continue '',$)') 200 ik=inkey(ik) if(ik.ne.13)goto 200 goto 1549 else goto 1565 endif endif c c...change amplitude c if(jj.eq.3)then write(*,'(1x// * 6x,'' New amplitude = '',$)') read(*,*,err=1549)temp if(temp.lt.0.d0)goto 1549 rifamp=temp goto 1549 endif c c...change phase c if(jj.eq.4)then write(*,'(1x// * 6x,'' New phase shift = '',$)') read(*,*,err=1549)temp if(temp.lt.-360.d0)goto 1549 if(temp.gt. 360.d0)goto 1549 rifph=temp goto 1549 endif c c...change shift c if(jj.eq.5)then write(*,'(1x// * 6x,'' New shift = '',$)') read(*,*,err=1549)temp rifsh=temp goto 1549 endif c c...subtract and exit c 1565 if(jj.eq.6.or.KK.eq.'i'.or.ncall.eq.0)then riffr=2.d0*pi*fifmhz riffr=riffr*tstep*1.d6 phase=rifph*2.d0*pi/360.d0 do 1555 j=1,nrep rif=rifamp*dcos(riffr*real(j)+phase)+rifsh idata(j)=idata(j)-rif 1555 continue ICHANG=1 if(ncall.eq.0)goto 1507 dummy=displaycursor($gcursoroff) goto 831 endif c dummy=displaycursor($gcursoroff) goto 831 C C C...display the interferogram help screen (with H) C 1600 dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) CALL clearscreen($GCLEARSCREEN) WRITE(*,1601) 1601 FORMAT(1x/'_____S I N G L E ', * ' I N T E R F E R O G R A M screen commands',20(1H_)/// * ' A/S - move lefthand cutoff cursor'/ * ' K/L - move righthand cutoff cursor'/ * ' caps on/off - fast/slow change in the above'/ * ' +/- - increase/decrease cursor step'// * ' B - background subtraction'/ * ' N - decay rate modification'/ * ' i/I - subtract remnant IF=20MHz (full AUTO/', * 'manual+AUTO)'/ * ' R - restore original interferogram'/ * ' U - ASCII dump of current interferogram'/ * ' O - subtraction of background interferogram'/ * ' up/dwn arrow - change coefficient for interferogram', * ' subtraction' //) WRITE(*,1602) 1602 FORMAT(17X, * 'Press ENTER to exit this HELP screen ',$) 1603 IK=INKEY(J) IF(IK.NE.13)GOTO 1603 CALL clearscreen($GCLEARSCREEN) GOTO 831 C C...subtraction of background interferogram (with O) C 1700 CALL clearscreen($GCLEARSCREEN) 1701 WRITE(*,1702) 1702 FORMAT(1x/'_____B A C K G R O U N D ', * 'I N T E R F E R O G R A M subtraction',14(1H_)/// * 5x,'The background interferogram is to be specified by a number', * ' defining '/ * 5x,'its position relative to the current interferogram'// * 5x,'..... ',$) read(*,'(i5)',err=1700)nbint if(nbint.eq.0)goto 1700 if(nmark+nbint.lt.0)goto 1700 if(nmark+nbint.gt.nscans)goto 1700 c nbint=ipoint(nmark+nbint) rmult=1.0 1704 indxi=intpt(nmark1,1)-1 indxj=intpt(nbint,1)-1 do 1703 j=1,nrep indxi=indxi+1 indxj=indxj+1 idata(j)=intall(indxi)-rmult*intall(indxj) 1703 continue c ichang=1 goto 831 C C...change scaling for background interferogram subtraction c (with up/down arrows) C 1750 if(ik.eq.-72)rmult=rmult*1.02 if(ik.eq.-80)rmult=rmult*0.980392 goto 1704 C c . . . . . . . . . . . . . . . . . . . . . . . . c c...display the help screen (with H) c 840 dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) CALL clearscreen($GCLEARSCREEN) WRITE(*,232) 232 FORMAT(1x/'_____F F T screen commands',53(1H_)/// * ' W/Z - change vertical scaling'/ * ' Q/E - change horizontal scaling'/ * ' A/S - shift spectrum left/right'/ * ' K/L/, - move marker left/right/centre'/ * ' caps on/off - fast/slow change in the above'// * ' I/ctrl I - show interferogram screen/interferogram info'/ * ' T - toggle between Power and Amplitude spectra'/ * ' N - change FFT zero filling parameter n'/ * ' P - show the FFT points'/ * ' R - rescale spectrum to initial conditions'/ * 10x,'U,Y - ASCII dump of current FFT: U is for standard FFT,'/ * 17x,'Y produces output scaled to mV range of interferogram'/ * 17x,'(use option Y only after having used option I)'/ * ' - exit back to previous screen'// * ' O - frequency of peak nearest the cursor'/ * ' 9 - use marker frequency as line frequency'/ * ' 0 - change height cutoff for peak measurement'/ * ' = - mean frequency of last two measured peaks'/ * ' - - intensity weighted mean frequency of last ', * 'two measured peaks'//) WRITE(*,106) 106 FORMAT(17X, * 'Press ENTER to exit this HELP screen ',$) 107 IK=INKEY(J) IF(IK.NE.13)GOTO 107 CALL clearscreen($GCLEARSCREEN) GOTO 699 c c...change FFT zero filling parameter n (with N) c 850 CALL clearscreen($GCLEARSCREEN) call FFTEXE <----- GOTO 851 c return end c C------------------------------------------------------------------------ c subroutine inpout(ioper,filnam,iexit,dirnam,ntsys) c c Routine to determine whether there are any spectral archives (extension c .FAR) in the directory and what is to be done with the spectra that are c found c c IOPER - operation code of what to do on return c = 0 read spectra according to LIST.DAT c = 1 read the archive with name in FILNAM c IEXIT = setting this to 1 on output requests termination by the M/P c DIRNAM = returns the name of the current directory on exit c USE DFLIB use globdat c logical*4 fsys parameter (maxarc=54,nivols=7) PARAMETER (ntextc=0, ntextb=7) character line*80,filnam*30,fnams(maxspe)*12,filarc*30 character cdummy*18,cdum*2,outa*12,outb*12,dirnam*50 integer*2 minfs(maxarc),maxfs(maxarc),intfre,iwk(maxspe), * iseen(maxspe),ndets(maxspe) real detvol(maxspe,2),volint(maxspe,nivols) real*8 wk(maxspe) real*8 fifmhz,fifkhz,fifmho,fifkho,fikhzs(maxspe) equivalence (cdum,intfre) common /scans/wk,wmult,iwk,filarc common /scans2/fikhzs,detvol,volint,iseen,nscans,ndets common /sfiles/fnams common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto C C C...determine the presence of FFT archives (with extension .FAR) C C this is done by means of a directory listing - note that default listing C format differs between W95/98 and NT/2000, and there are also small C differences between language versions. C C The DIR keyword -N enforces in W2000 the old 8.3 form of listing, but C it is an illegal option in W95/98 C ntsys=0 fsys=systemqq('ver>far.lst') open (2,file='far.lst',status='old') read(2,'(a)',end=9)line read(2,'(a)',end=9)line if(line(1:1).eq.'M')ntsys=1 close(2) c iexit=0 write(*,'(1x/)') if(ntsys.eq.1)then fsys=systemqq('dir *.far/On/-N>far.lst') else fsys=systemqq('dir *.far/On>far.lst') endif if(fsys.neqv..TRUE.)then write(*,'(1x/a//)')' ***** ERROR: Cannot do dir *.far>far.lst' stop endif open(2,file='far.lst',status='old') c c...skip top four lines and determine directory name from the fifth line c do 150 i=1,4 read(2,'(a)',end=9)line 150 continue i=len_trim(line) if(line(2:9).eq.'Katalog ')then dirnam=line(10:i) endif if(line(2:9).eq.'Katalog:')then dirnam=line(11:i) endif if(line(2:14).eq.'Directory of ')dirnam=line(15:i) read(2,'(a)',end=9)line c narch=0 7 read(2,'(a)',end=9)line if(line(10:12).ne.'FAR'.and.line(10:12).ne.'far')goto 7 do 8 i=9,1,-1 if(line(i:i).ne.' ')goto 10 8 continue 10 narch=narch+1 c if(narch.eq.maxarc)then dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) write(*,'(1x//'' ***** Too many archives - the program is '', * ''only dimensioned up to'',i5//)')maxarc dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) stop endif fnams(narch)=line(1:i)//'.'//line(10:12) goto 7 c 9 close(2) fsys=systemqq('del far.lst') c c...options c 19 if(narch.eq.0)then write(*,100) 100 format(1x// * ' ----> The directory '$) DUMmy2=settextcolor(int2(12)) write(*,'(a)')dirnam(1:len_trim(dirnam)) DUMmy2=settextcolor(int2(ntextc)) write(*,'('' contains NO FFT archive files - '', * ''attempting to use the LIST.DAT file'')') 101 ioper=0 return else i=len_trim(dirnam) DUMmy2=settextcolor(int2(12)) write(*,'(1x,a,$)')dirnam(1:i) DUMmy2=settextcolor(int2(ntextc)) write(*,15)narch 15 format(' contains',i3,' FFT archive(s):'/40(1H )) c c...read frequency limits of archives c do 200 i=1,narch iwk(i)=i open(2,file=fnams(i),form='binary',status='old',err=200) read(2)cdummy read(2)cdum if(cdum.ne.'--')then minfs(i)=intfre wk(i)=intfre else wk(i)=0.0d0 endif read(2)cdum if(cdum.ne.'--')maxfs(i)=intfre close(2) 200 continue c c...sort according to frequency (-ve NSCANS is used to avoid c further reordering built into this version of SORTH) c nscans=-narch if(narch.gt.1)call sorth <----- nscans=0 c c...write list of archives c c...single column c if(narch.lt.10)then do 12 j=1,narch i=iwk(j) if(minfs(i).eq.0.and.maxfs(i).eq.0)then write(outa,'(a12)')' ' else write(outa,'(i6,''-'',i5)')minfs(i),maxfs(i) endif write(*,17)j,fnams(i),outa 12 continue 17 format(12x,i2,' = ',a,a) endif c c...two column c if(narch.le.36.and.narch.ge.10)then j=narch/2 c do 112 k=1,j line by line i=iwk(k) ipj=iwk(k+j) if(minfs(i).eq.0.and.maxfs(i).eq.0)then write(outa,'(a12)')' ' else write(outa,'(i6,''-'',i5)')minfs(i),maxfs(i) endif if(minfs(ipj).eq.0.and.maxfs(ipj).eq.0)then write(outb,'(a12)')' ' else write(outb,'(i6,''-'',i5)')minfs(ipj),maxfs(ipj) endif write(*,117)k,fnams(i),outa,k+j,fnams(ipj),outb 112 continue c if(2*(narch/2).ne.narch)then last odd file i=iwk(narch) if(minfs(i).eq.0.and.maxfs(i).eq.0)then write(outa,'(a12)')' ' else write(outa,'(i6,''-'',i5)')minfs(i),maxfs(i) endif write(*,118)narch,fnams(i),outa endif c 117 format(4x,i2,' = ',a,a, 12x,i2,' = ',a,a) 118 format(45x, i2,' = ',a,a) endif c c...three column c if(narch.gt.36.and.narch.le.maxarc)then j=narch/3 c do 120 k=1,j i=iwk(k) ipj=iwk(k+j) iipj=iwk(k+2*j) write(*,121)k,fnams(i),k+j,fnams(ipj),k+2*j,fnams(iipj) 120 continue c if(narch-j*3.eq.1)then i=iwk(narch) write(*,122)narch,fnams(i) endif c if(narch-j*3.eq.2)then i=iwk(narch-1) write(*,122)narch-1,fnams(i) i=iwk(narch) write(*,122)narch,fnams(i) endif c 121 format(5x, i2,' = ',a,2(10x,i2,' = ',a)) 122 format(59x,i2,' = ',a) endif c c write(*,18) 18 format(1x/ * ' OPTIONS: 0 = read spectra according to summary data file'/ * ' n = read archive n (any -ve value for exit)'// * 25x,'..... '$) read(*,'(i2)',err=20)ioper if(ioper.lt.0)then iexit=1 return endif if(ioper.lt.0.or.ioper.gt.narch)goto 20 if(ioper.gt.0)then filnam=fnams(iwk(ioper)) ioper=1 else ioper=0 endif return endif c 20 goto 19 c return end c C_____________________________________________________________________________ c subroutine inpspe(nfil,ntsys) c c Routine to set up a list of files available for input as spectra c c USE DFLIB c logical*4 fsys parameter (maxspe=500) PARAMETER (ntextc=0, ntextb=7) character line*80,fnams(maxspe)*12 common /sfiles/fnams c c c...save the current directory to file c write(*,'(1x/)') if(ntsys.eq.1)then fsys=systemqq('dir /Od/-N>spec.lst') else fsys=systemqq('dir /Od>spec.lst') endif if(fsys.neqv..TRUE.)then write(*,'(1x/a//)')' ***** ERROR: Cannot do dir >spec.lst' stop endif c c...go through the directory file and identify potential spectral files c open(2,file='spec.lst',status='old') c do 150 i=1,5 read(2,'(a)',end=9)line 150 continue c nfil=0 c 7 read(2,'(a)',end=9)line c 21 if(line(1:1).eq.' '.or.line(1:1).eq.'.')goto 7 if(line(10:10).eq.'~')goto 7 if(line(10:12).eq.'FAR'.or.line(10:12).eq.'far')goto 7 if(line(10:12).eq.'FOR'.or.line(10:12).eq.'for')goto 7 if(line(10:12).eq.'EXE'.or.line(10:12).eq.'exe')goto 7 if(line(10:12).eq.'OBJ'.or.line(10:12).eq.'obj')goto 7 if(line(10:12).eq.'DAT'.or.line(10:12).eq.'dat')goto 7 if(line(10:12).eq.'LST'.or.line(10:12).eq.'lst')goto 7 if(line(14:14).eq.'<'.or.line(16:16).eq.'<')goto 7 if(line(10:12).eq.'OUT'.or.line(10:12).eq.'out')goto 7 if(line(10:12).eq.'OUT'.or.line(10:12).eq.'out')goto 7 if(line(10:12).eq.'ARJ'.or.line(10:12).eq.'arj')goto 7 if(line(10:12).eq.'ASR'.or.line(10:12).eq.'asr')goto 7 if(line(10:12).eq.'INP'.or.line(10:12).eq.'inp')goto 7 if(line(10:12).eq.'ASF'.or.line(10:12).eq.'asf')goto 7 if(line(10:12).eq.'PAR'.or.line(10:12).eq.'par')goto 7 if(line(10:12).eq.'LIN'.or.line(10:12).eq.'lin')goto 7 if(line(10:12).eq.'VAR'.or.line(10:12).eq.'var')goto 7 if(line(10:12).eq.'BIN'.or.line(10:12).eq.'bin')goto 7 if(line(10:12).eq.'INT'.or.line(10:12).eq.'int')goto 7 if(line(10:12).eq.'FIT'.or.line(10:12).eq.'fit')goto 7 if(line(10:12).eq.'CAT'.or.line(10:12).eq.'cat')goto 7 if(line(10:12).eq.'PMI'.or.line(10:12).eq.'pmi')goto 7 if(line(10:12).eq.'COR'.or.line(10:12).eq.'cor')goto 7 if(line(10:12).eq.'STF'.or.line(10:12).eq.'stf')goto 7 if(line(10:12).eq.'GLE'.or.line(10:12).eq.'gle')goto 7 if(line(10:12).eq.'DOC'.or.line(10:12).eq.'doc')goto 7 if(line(10:12).eq.'BAK'.or.line(10:12).eq.'bak')goto 7 if(line(1:6).eq.'STATUS'.and.line(10:11).eq.'ME')goto 7 c do 8 i=9,1,-1 if(line(i:i).ne.' ')goto 10 8 continue 10 nfil=nfil+1 fnams(nfil)=line(1:i)//'.'//line(10:12) if(nfil.eq.maxspe)then dummy4=setbkcolor(12) DUMmy2=settextcolor(int2(15)) write(*,'(1x//'' ***** The current limit of'',i3, * '' on the number of files has been reached''/ * '' ***** NO MORE WILL BE READ IN''/)') * maxspe dummy4=setbkcolor(ntextb) DUMmy2=settextcolor(int2(ntextc)) goto 9 endif goto 7 9 close(2) fsys=systemqq('del spec.lst') c return end c C_____________________________________________________________________________ c subroutine FFTEXE c C Intermediate routine which prepares input for the FFT calculation c on the recorded interferogram. C C The amount of zero-filling is determined by the variable MULZER - this c can be set externally to this routine (but given a negative sign), in which c case no screen output is made by this routine. C use globdat PARAMETER (Nmaxpt=131072) c real*4 data(nmaxpt),w1(2*nmaxpt),w2(nmaxpt) real*4 p(nmaxpt) integer*4 idata(maxpts) real*8 fstep logical ovrlap common /points/data,npts common /work/w1 common /work1/w2 common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /specf/p,fstep,ncut,MULZER,NCALL C c c...Transfer points to array DATA passed to FFT routines: c minimum number of points is 256 or higher powers of 2 so fill-up c appropriately c NCALL=NCALL+1 NOWRIT=0 i=0 do 170 j=nskips+1,nrep-nskipe i=i+1 data(i)=idata(j) 170 continue npts=i ifromi=i IF(NPTS.LE.131072)NMAX=131072 IF(NPTS.LE. 65536)NMAX= 65536 IF(NPTS.LE. 32768)NMAX= 32768 IF(NPTS.LE. 16384)NMAX= 16384 IF(NPTS.LE. 8192)NMAX= 8192 if(npts.le. 4096)nmax= 4096 if(npts.le. 2048)nmax= 2048 if(npts.le. 1024)nmax= 1024 if(npts.le. 512)nmax= 512 if(npts.le. 256)nmax= 256 if(npts.ne.nmax)then do 52 i=npts+1,nmax data(i)=data(npts) 52 continue npts=nmax endif c c...FFT options: amount of zero-filling: note that -ve MULZER can be C generated by routine LOOKIN for recalculation of FFT with previous c value of MULZER, ie. when cutoffs for rejected points have been modified C C -on first time into this routine (NCALL=1) default value of MULZER is c assigned c IF(MULZER.GE.0)THEN write(*,56)mulzer,ifromi,npts,npts*2**mulzer 56 FORMAT(1x/'_____C H A N G E F F T Z E R O F I L L I N G', * 30(1H_)/// * 5x,'The last FFT was performed using filling parameter n=',i1, * ' and:'// * i7,' data points from the interferogram, first filled-up to'/ * i7,' points (= Npts), and then to a total of'/ * i7,' points by replicating the last data point.'/) 31 if(NCALL.ne.1)then 33 write(*,30) 30 format(1x/5x,'Specify new value of the filling parameter n ', * 'defined by'/ * 5x,'Total_Npts = 2**n * Npts (n=0,1,2,3.... )'// * 45x,' n = ',$) read(*,*,err=31)mulzer if(npts*2**mulzer.lt.npts)goto 33 else mulzer=4 555 if(npts*2**mulzer.gt.nmaxpt)then MULZER=MULZER-1 GOTO 555 ENDIF ENDIF C if(npts*2**mulzer.gt.nmaxpt)then mulzer=npts*2**mulzer write(*,775)mulzer,nmaxpt 775 format(1x/' ERROR:',i8,' points needed for FFT but only ', * i6,' dimensioned'// * ' -----> Specify smaller n'/) goto 31 endif ELSE mulzer=-mulzer 556 if(npts*2**mulzer.gt.nmaxpt)then MULZER=MULZER-1 GOTO 556 ENDIF NOWRIT=1 ENDIF C SHIFT=DATA(NPTS) do 32 i=1,npts*2**mulzer IF(I.LE.NPTS)THEN DATA(I)=DATA(I)-SHIFT ELSE DATA(I)=0.0 ENDIF 32 continue npts=npts*2**mulzer C k=1 m=npts/(k+1) IF(NOWRIT.NE.1)write(*,57)npts 57 format(1x/i7,' point record will be used for FFT') c c...FFT c samfre=1./tstep fnyq=samfre/2. fstep=samfre/npts ncut=fnyq/fstep fstep=fstep/1000. IF(NOWRIT.NE.1) * write(*,25)tstep*1.E+6,samfre*1.E-6,fnyq*1.E-6,fstep,ncut 25 format(1x/' time step = ',f15.10,' microsec.'/ * ' sampling frequency = ',f15.10,' MHz'/ * ' Nyquist frequency = ',f15.10,' MHz'/ * ' frequency step = ',f15.10,' kHz'/ * ' points in Nyq. interval = ',i15/) c ovrlap=.true. call spctrm(m,k,ovrlap) <----- C C...Postscale points in power spectrum so that their intensities for c various amounts of zero filling are unified and equivalent to those c for n=4 c ymult=mulzer-4 ymult=16.d0**nint(ymult) if(nint(ymult).ne.1)then do 570 i=1,ncut p(i)=p(i)*ymult 570 continue endif c return end C C------------------------------------------------------------------------ c SUBROUTINE spctrm(m,k,ovrlap) c c Power spectrum estimation using routine 'four1' c c p = on output contains the input data's power (mean square amplitude) at c frequency (j-1)/(2*m) cycles per gridpoint, for j=1,2....,m c m = number of data points in segment c k = number of segments (each with 2m data points) c ovrlap=.false. segments do not overlap, 4*m*k data points c ovrlap=.true. segments overlap, (2k+1)*m data points c data = time domain data points c parameter (nmax=131072) INTEGER*4 k,m c REAL p(m),w1(2*nmax),w2(nmax) REAL*4 p(nmax),w1(2*nmax),w2(nmax) LOGICAL ovrlap real*4 data(nmax) real*8 fstep common /points/data,npts common /work/w1 common /work1/w2 common /specf/p,fstep,ncut,NFFT,NCALL c INTEGER*4 j,j2,joff,joffn,kk,m4,m43,m44,mm REAL*4 den,facm,facp,sumw,w,window c window(j)=(1.-abs(((j-1)-facm)*facp)) Bartlett c window(j)=1. Square c window(j)=(1.-(((j-1)-facm)*facp)**2 Welch nread=0 mm=m+m m4=mm+mm m44=m4+4 m43=m4+3 den=0. facm=m facp=1./m c c...accumulate the squared sum of the weights c sumw=0. do 11 j=1,mm sumw=sumw+window(j)**2 11 continue c c...initialize the spectrum to zero c do 12 j=1,m p(j)=0. 12 continue c c...initialize the 'save' half-buffer - this is a modifcation to use c the data from common block /points/. The values are read in c successively and NREAD is the total number of data points used. c If more points are required then are in the data then the last point c is repeated c if(ovrlap)then do 21 j=1,m nread=nread+1 if(nread.gt.npts)then w2(j)=data(npts) else w2(j)=data(nread) endif 21 continue endif c c...Loop over data set segments in groups of two. Get two complete c segments into workspace. c do 18 kk=1,k do 15 joff=-1,0,1 if(ovrlap)then do 13 j=1,m w1(joff+j+j)=w2(j) 13 continue do 22 j=1,m nread=nread+1 if(nread.gt.npts)then w2(j)=0. else w2(j)=data(nread) endif 22 continue joffn=joff+mm do 14 j=1,m w1(joffn+j+j)=w2(j) 14 continue else do 23 j=joff+2,m4,2 nread=nread+1 if(nread.gt.npts)then w1(j)=0. else w1(j)=data(nread) endif 23 continue endif 15 continue c c...Apply the window to the data c do 16 j=1,mm j2=j+j w=window(j) w1(j2)=w1(j2)*w w1(j2-1)=w1(j2-1)*w 16 continue c c...Fourier transform the windowed data c call four1(mm,1) <----- c c...Sum results into previous segments c p(1)=p(1)+w1(1)**2+w1(2)**2 do 17 j=2,m j2=j+j p(j)=p(j)+w1(j2)**2+w1(j2-1)**2 * +w1(m44-j2)**2+w1(m43-j2)**2 17 continue den=den+sumw 18 continue c c...Correct normalization and normalize the output c den=m4*den do 19 j=1,m p(j)=p(j)/den 19 continue c write(*,25)nread 25 format(i10,' points used in FFT') c return end c c---------------------------------------------------------------------------- c SUBROUTINE four1(nn,isign) parameter (nmax=131072) INTEGER*4 isign,nn c REAL data(2*nn) real*4 data(2*nmax) common /work/data c c Routine replaces data(1:2*nn) by its discrete Fourier transform, if isign c is input as -1; or replaces data(1:2*nn) by nn times its inverse discrete c Fourier transform, if isign is input as -1. c data is a complex array of length nn, or equivalently, a real array of c length 2*nn. c nn MUST be an integer power of 2 (this is not checked for!) c INTEGER*4 i,istep,j,m,mmax,n REAL*4 tempi,tempr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp n=2*nn j=1 do 11 i=1,n,2 if(j.gt.i)then tempr=data(j) tempi=data(j+1) data(j)=data(i) data(j+1)=data(i+1) data(i)=tempr data(i+1)=tempi endif m=n/2 1 if((m.ge.2).and.(j.gt.m))then j=j-m m=m/2 goto 1 endif j=j+m 11 continue c mmax=2 2 if(n.gt.mmax)then istep=2*mmax theta=6.28318530717959d0/(isign*mmax) wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0 wi=0.d0 do 13 m=1,mmax,2 do 12 i=m,n,istep j=i+mmax tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1) tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j) data(j)=data(i)-tempr data(j+1)=data(i+1)-tempi data(i)=data(i)+tempr data(i+1)=data(i+1)+tempi 12 continue wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi 13 continue mmax=istep goto 2 endif c return end c C------------------------------------------------------------------------ c SUBROUTINE PF(fmark,smax,XPEAK,ERRORX,htcut,xbot,fhalm,fhalp) C C Position of line maximum is established using the bisection method: C midpoints of line contour are determined at preselected number of sections c NSECT from peak maximum down to a selected fraction of peak height, HTCUT. c C Straight line fit to such midpoints gives, for y equal to line maximum, c the line frequency, with some account for possible line asymmetry. C C FMARK - frequency of the marker, which is assumed to have been set near C peak maximum c SMAX - on exit the intensity at the maximum C XPEAK - on exit the required central fitted peak position (requires C addition of FCENT) C ERRORX - on exit error on the fitted peak position c HTCUT - proportion of line height to which profile division is taken c XBOT - value of X for Y=0.4Ymax for drawing fitted line c FHALM - offset for negative FWHH point c FHALP - offset for positive FWHH point C PARAMETER (Nmaxpt=131072,nsect=20) c real p(nmaxpt) real*8 fstep,fmark,x(nsect),y(nsect),fr(nsect,2),xpeak REAL*8 SUMX,SUMY,SUMXY,SUMX2,SUMY2,CXX,CXY,CYY,RN,A0,A1, * xoffs,yoffs c common /specf/p,fstep,npts,NFFT,NCALL common /peak/x,y C c...determine initial value of peak maximum and its position c n=(fmark/fstep)+1. SMAX=P(N) 1 IF(N.LE.1.OR.N.GE.NPTS-1)THEN XPEAK=(N-1)*FSTEP ERRORX=(NPTS-1)*FSTEP RETURN ENDIF IF(P(N+1).GT.SMAX)THEN SMAX=P(N+1) N=N+1 GOTO 1 ENDIF IF(P(N-1).GT.SMAX)THEN SMAX=P(N-1) N=N-1 GOTO 1 ENDIF NMAX=N C C...Determine frequencies of points on sections through line contour: c linear interpolation used C do 2 ns=1,nsect ysect=smax-ns*smax*htcut/real(nsect) y(ns)=ysect do 3 n=nmax,npts if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fr(ns,2)=(rn-1.0)*fstep goto 2 endif 3 continue 2 continue c do 5 ns=1,nsect ysect=y(ns) do 4 n=nmax,2,-1 if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fr(ns,1)=(rn-1.0)*fstep goto 5 endif 4 continue 5 continue C do 6 n=1,nsect x(n)=0.5*(fr(n,1)+fr(n,2)) 6 continue c C...Straight line fit: since the line is almost vertical (gradient very c large) this was found to lead to numerical instabilities and for c this reason axes are reversed for least squares. For further increased c numerical stability SMAX is subtracted from Y and X(1) from X ie. the c equation of fit is: c c (x-x1) = a0 + a1 (y-smax) C SUMX=0.D0 SUMY=0.D0 SUMXY=0.D0 SUMX2=0.D0 SUMY2=0.D0 XOFFS=X(1) YOFFS=SMAX DO 7 I=1,NSECT SUMy=SUMy+(X(I)-xoffs) SUMx=SUMx+(Y(I)-yoffs) SUMXY=SUMXY+(X(I)-xoffs)*(Y(I)-yoffs) SUMy2=SUMy2+(X(I)-xoffs)**2 SUMx2=SUMx2+(Y(I)-yoffs)**2 7 CONTINUE C C...coefficients RN=NSECT CXX=SUMX2-SUMX*SUMX/RN CXY=SUMXY-SUMX*SUMY/RN A1=CXY/CXX IF(A1.EQ.0.D0)THEN ERRORX=0.D0 XPEAK=0.D0 RETURN ENDIF A0=(SUMY-A1*SUMX)/RN c c...peak frequency xpeak=a0+x(1) c C...coordinates of 0.4Ymax point xbot=a0-0.6*smax*a1+x(1) C C...error CYY=SUMY2-SUMY*SUMY/RN ERA1S=((CYY/CXX)-(CXY/CXX)**2)/(RN-2.D0) ERA0S=SUMX2*ERA1S/RN ERRORX=dsqrt(dble(era0s)) C C...Find X values at FWHH C ysect=0.5*smax do 13 n=nmax,npts if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fhalp=(rn-1.0)*fstep goto 12 endif 13 continue fhalp=(npts-1)*fstep c 12 do 14 n=nmax,2,-1 if(p(n).le.ysect)then rn=real(n-1)+(ysect-p(n-1))/(p(n)-p(n-1)) fhalm=(rn-1.0)*fstep goto 15 endif 14 continue fhalm=fstep 15 continue c RETURN END C C------------------------------------------------------------------------ C SUBROUTINE SORTH c c This routine is based on the SORT2 'heapsort' routine from Numerical c Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N) C in ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities. c use globdat c parameter (nivols=7) character fnams(maxspe)*12,ftemp(maxspe)*12,filarc*30 real detvol(maxspe,2),volint(maxspe,nivols) real*8 fifmhz,fifkhz,fifmho,fifkho,fikhzs(maxspe) integer*2 ipt(maxspe),iseen(maxspe),ndets(maxspe) real*8 wk(maxspe) common /scans/wk,wmult,ipt,filarc common /scans2/fikhzs,detvol,volint,iseen,nscans,ndets common /sfiles/fnams common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto c INTEGER*2 IIPT,L,N,NSTART,I,J,IR REAL*8 WWK real rtemp(maxspe) equivalence (rtemp(1),ftemp(1)) C nstart=1 n=iabs(nscans) c L=N/2+1 IR=N 10 CONTINUE IF(L.GT.NSTART)THEN L=L-1 WWK=WK(L) IIPT=IPT(L) ELSE WWK=WK(IR) IIPT=IPT(IR) WK(IR)=WK(1) IPT(IR)=IPT(1) IR=IR-1 IF(IR.EQ.NSTART)THEN WK(1)=WWK IPT(1)=IIPT GOTO 100 ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(WK(J).LT.WK(J+1))J=J+1 ENDIF IF(WWK.LT.WK(J))THEN WK(I)=WK(J) IPT(I)=IPT(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF WK(I)=WWK IPT(I)=IIPT GO TO 10 c c...reorder FNAMS, DETVOL and VOLINT according to the order in IPT c (this is not done if NSCANS has previously been made negative) c 100 if(nscans.lt.0)return do 101 l=1,2 do 102 i=1,nscans rtemp(i)=detvol(i,l) 102 continue do 101 i=1,nscans j=ipt(i) detvol(i,l)=rtemp(j) 101 continue c do 201 l=1,nivols do 202 i=1,nscans rtemp(i)=volint(i,l) 202 continue do 201 i=1,nscans j=ipt(i) volint(i,l)=rtemp(j) 201 continue c do 302 i=1,nscans ftemp(i)=fnams(i) 302 continue do 301 i=1,nscans j=ipt(i) fnams(i)=ftemp(j) 301 continue c c RETURN END C C_____________________________________________________________________________ c subroutine baksub(nspt) c use globdat PARAMETER (Nmaxpt=131072,maxsmo=1001) c integer*4 idata(maxpts),ioldat(maxpts),itemp(maxpts) real*4 spol(maxsmo) common /scans1/idata,vstep,tstep,nrep,nskips,nskipe,naver common /smooth/ioldat,itemp,spol c C...Subtraction of background by triple smoothing of interferogram with c least squares smoothing interval c C For smoothing interval of length 2m+1 the elements of the smoothing c (cubic) polynomial are given by: C C 3(3m**2 + 3m -1 - 5s**2) C c(s) = ------------------------ C (2m+1) (2m-1) (2m+3) C C where s runs from -m to +m (T.H.Edwards and P.D.Wilson, Applied C Spectroscopy 28,541-545(1974)) C C c...set up coefficients in smoothing polynomial M=NSPT/2 T1=3.D0/((2*M+1)*(2.D0*M-1.D0)*(2*M+3)) T2=3*M*M+3.D0*M-1.D0 DO 1103 j=1,NSPT IS=j-M-1 SPOL(j)=T1*(T2-5*IS*IS) 1103 CONTINUE C c...Smooth three times an ITEMP copy of data in IDATA into IDATA c ISTRT=M+1 IFIN=nrep-M DO 1104 k=1,3 do 1105 j=1,nrep itemp(j)=idata(j) 1105 continue DO 543 I=1,nrep SUM=0. DO 544 J=1,NSPT IS=J-M-1 II=I+IS IF(II.LT.1)II=iabs(II)+1 IF(II.GT.nrep)II=nrep-(II-nrep-1) SUM=SUM+itemp(II)*SPOL(J) 544 CONTINUE Idata(I)=sum 543 CONTINUE DO 545 I=1,nrep ITEMP(I)=IDATA(I) 545 CONTINUE 1104 continue C c...subtract baseline in IDATA from original interferogram in IOLDAT c do 1106 j=1,nrep idata(j)=ioldat(j)-idata(j) 1106 continue c return end C C_____________________________________________________________________________ c integer*2 function INKEY(N2) c c By L Pszczolkowski: c c...This emulates for MSF PS1.0 the INKEY function of Z.Czumaj which c in turn emulated for MSF5.0 the INKEY function from IIUWGRAF graphics c library for the Hercules card c c The function GETCHARQQ returns the ASCII character if the corresponding c key was pressed. If function or direction key was pressed then 0 c or hex E0 is returned and another call to GETCHARQQ is required to c get the extended code of the character c USE DFLIB c INTEGER*2 IK CHARACTER*1 KK c KK=GETCHARQQ() IK=ICHAR(KK) IF(IK.EQ.0 .OR. IK.EQ.224 ) THEN KK=GETCHARQQ() IK=-ICHAR(KK) ENDIF n2=ik INKEY=IK END C_____________________________________________________________________________ C C Non-linear fit C_____________________________________________________________________________ C C SUBROUTINE MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT, 1 COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA) C C............................................................................. C C Levenberg-Marquardt method, attempting to reduce the value CHISQ of a fit C between a set of NDATA points X(I),Y(I) with individual standard devia- C tions SIG(I), and a nonlinear function dependent on MA coefficients A. C C C On the first call of MRQMIN provide an initial guess for the parameters C A, and set ALAMDA<0 for initialization (which then sets ALAMDA=0.001). C If a step succeeds CHISQ becomes smaller and ALAMDA decreases C by a factor of 10. C If a step fails ALAMDA grows by a factor of 10. You must call this C routine repeatedly until convergence is achieved. Then make one final C call with ALAMDA=0., so that COVAR(I,J) returns the covariance matrix, C and ALPHA(I,J) the curvature matrix. C C NOTE: This is double precision version of the routine set from Num.Rec. C and there is also an extra CALL SETUP statement in MRQCOF to C work with STRFIT - these and other modifications identified C with ...ZK C............................................................................. C C ROUTINES REQUIRED: C C MRQCOF, COVSRT, GAUSSJ - called by MRQMIN C FUNCS(user supplied) - called by MRQCOF C C C X,Y,SIG - vectors of length NDATA containing X,Y pairs and deviations C on the Y values C A - vector of length MA containing the constants to be used for C evaluations of Y.calc (MFIT of which are to be fitted) C LISTA - vector of length MFIT carrying pointers to the constants in C A which are to be fitted C COVAR,ALPHA - square arrays with physical dimension NCA (which is to be C greater or equal to MFIT) used as workspace. Routine C COVSRT will return the covariance matrix in COVAR. C CHISQ - minimisation parameter chi-square C FUNCS - user supplied subroutine FUNCS(X,A,YFIT,DYDA,MA) that C evaluates the calculated value of Y: YFIT at X, and the C derivatives of Y with respect to the fitting parameters A C at X. C ALAMDA - parameter to monitor the prograss of fit, to be set to any C -ve value for initialisation. This will then be set to 0.001 C and on following cycles decreased by a factor of 10 is chi- C square decreases and increased by a factor of 10 if CHISQ C increases. ALAMDA=0 returns the covariance matrix COVAR C and the curvature (Hessian) matrix ALPHA C............................................................................. C C...MMAX defines the largest number of parameters C IMPLICIT REAL*8 (A-H,O-Z) ...ZK EXTERNAL FUNCS ...ZK PARAMETER (maxpar=10, MMAX=maxpar) ...ZK COMMON /COV/COVTMP(maxpar,maxpar) ...ZK DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MAXPAR),LISTA(MAXPAR), ...ZKmod 1 COVAR(NCA,NCA),ALPHA(NCA,NCA),ATRY(MMAX),BETA(MMAX),DA(MMAX) c c...initialization c IF(ALAMDA.LT.0.d0)THEN KK=MFIT+1 c c...does LISTA contain a proper permutation of the coefficients? DO 12 J=1,MA IHIT=0 DO 11 K=1,MFIT IF(LISTA(K).EQ.J)IHIT=IHIT+1 11 CONTINUE IF (IHIT.EQ.0) THEN LISTA(KK)=J KK=KK+1 ELSE IF (IHIT.GT.1) THEN PAUSE 'Improper permutation in LISTA' ENDIF 12 CONTINUE IF (KK.NE.(MA+1)) PAUSE 'Improper permutation in LISTA' ALAMDA=0.001d0 CALL MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NCA,CHISQ, 1 FUNCS) <----- OCHISQ=CHISQ DO 13 J=1,MA ATRY(J)=A(J) 13 CONTINUE ENDIF c c...alter linearized fitting matrix, by augmenting diagonal elements c DO 15 J=1,MFIT DO 14 K=1,MFIT COVAR(J,K)=ALPHA(J,K) 14 CONTINUE COVAR(J,J)=ALPHA(J,J)*(1.+ALAMDA) DA(J)=BETA(J) 15 CONTINUE c c...matrix solution c CALL GAUSSJ(COVAR,MFIT,NCA,DA,1,1) <----- c c...once converged evaluate covariance matrix with ALAMDA=0 c (but if ALAMDA not zero then transfer covariance matrix to ...ZK c array COVTMP for passing over to calling program for evaluation ...ZK c of errors for each cycle) ...ZK c IF(ALAMDA.NE.0.d0)THEN ...ZK DO 20 IZA=1,MA ...ZK DO 20 IZB=1,MA ...ZK COVTMP(IZA,IZB)=COVAR(IZA,IZB) ...ZK 20 CONTINUE ...ZK ENDIF IF(ALAMDA.EQ.0.d0)THEN CALL COVSRT(COVAR,NCA,MA,LISTA,MFIT) <----- RETURN ENDIF c c...did the trial succeed? c DO 16 J=1,MFIT ATRY(LISTA(J))=ATRY(LISTA(J))+DA(J) 16 CONTINUE CALL MRQCOF(X,Y,SIG,NDATA,ATRY,MA,LISTA,MFIT,COVAR,DA,NCA,CHISQ, 1 FUNCS) <----- c c...success, accept the new solution c IF(CHISQ.LT.OCHISQ)THEN ALAMDA=0.1d0*ALAMDA OCHISQ=CHISQ DO 18 J=1,MFIT DO 17 K=1,MFIT ALPHA(J,K)=COVAR(J,K) 17 CONTINUE BETA(J)=DA(J) A(LISTA(J))=ATRY(LISTA(J)) 18 CONTINUE c c...failure, increase ALAMDA and return c ELSE ALAMDA=10.d0*ALAMDA C C the following has been commented out to avoid termination of iterations C CHISQ=OCHISQ ...ZK ENDIF c RETURN END c c_____________________________________________________________________________ c SUBROUTINE MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NALP, 1 CHISQ,FUNCS) C C Routine used by MRQMIN to evaluate the linearized fitting matrix ALPHA, C and vector BETA (equation 14.4.8 Num.Rec.) C C IMPLICIT REAL*8 (A-H,O-Z) ...ZK PARAMETER (maxpar=10, MMAX=maxpar) ...ZK DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),ALPHA(NALP,NALP),BETA(MA), * DYDA(MMAX),LISTA(MAXPAR),A(MAXPAR) ...ZKmod c c...initialize (symmetric) ALPHA,BETA c DO 12 J=1,MFIT DO 11 K=1,J ALPHA(J,K)=0.d0 11 CONTINUE BETA(J)=0.d0 12 CONTINUE CHISQ=0. c c...summation loop over all data C DO 15 I=1,NDATA CALL FUNCS(X(I),A,YMOD,DYDA,MA) <--FUNCS SIG2I=1.d0/(SIG(I)*SIG(I)) DY=Y(I)-YMOD DO 14 J=1,MFIT WT=DYDA(LISTA(J))*SIG2I DO 13 K=1,J ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K)) 13 CONTINUE BETA(J)=BETA(J)+DY*WT 14 CONTINUE c c...find CHISQ CHISQ=CHISQ+DY*DY*SIG2I 15 CONTINUE c c...fill in the symmetric side c DO 17 J=2,MFIT DO 16 K=1,J-1 ALPHA(K,J)=ALPHA(J,K) 16 CONTINUE 17 CONTINUE c RETURN END c c_____________________________________________________________________________ c SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) C C Routine for solving a set of linear equations by the Gauss-Jordan C elimination with full pivoting. C C A is an input matrix of N by N elements, stored in an array of physical C dimensions NP by NP. C B is an input matrix of N by M containing the M right-hand side vectors, C stored in an array of physical dimensions NP by MP. C On output, A is replaced by its matrix inverse, and B is replaced by the C corresponding set of solution vectors C C IMPLICIT REAL*8 (A-H,O-Z) ...ZK PARAMETER (NMAX=50) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) c DO 11 J=1,N IPIV(J)=0 11 CONTINUE c c...this is the main loop over the columns to be reduced c DO 22 I=1,N BIG=0.d0 c c...this is the outer loop of the search for the pivot element c DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'Singular matrix' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 c c...we now have the pivot element, so we interchange rows, if needed to put c the pivot on the diagonal. The columns are not physically interchanged, c only relabeled: c IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.d0) PAUSE 'Singular matrix.' PIVINV=1.d0/A(ICOL,ICOL) A(ICOL,ICOL)=1.d0 DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0.d0 DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE c c...this is the end of the main loop over columns of the reduction. c It only remains to unscramble the solution in view of the column c interchanges. c DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE c RETURN END c c_____________________________________________________________________________ c SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT) C C Given the covariance matrix COVAR of a fit for MFIT of MA total parameters C and their ordering LISTA(I), repack the covariance matrix to the true C order of the parameters. C Elements associated with fixed parameters will be zero. NCVM is the C physical dimension of COVAR C C IMPLICIT REAL*8 (A-H,O-Z) ...ZK PARAMETER (maxpar=10) DIMENSION COVAR(NCVM,NCVM),LISTA(MAXPAR) ...ZKmod c c...zero all elements below diagonal c DO 12 J=1,MA-1 DO 11 I=J+1,MA COVAR(I,J)=0.d0 11 CONTINUE 12 CONTINUE c c...repack off-diagonal elements of fit into correct locations below diagonal c DO 14 I=1,MFIT-1 DO 13 J=I+1,MFIT IF(LISTA(J).GT.LISTA(I)) THEN COVAR(LISTA(J),LISTA(I))=COVAR(I,J) ELSE COVAR(LISTA(I),LISTA(J))=COVAR(I,J) ENDIF 13 CONTINUE 14 CONTINUE c c...temporarily store original diagonal elements in top row, and zero the c diagonal c SWAP=COVAR(1,1) DO 15 J=1,MA COVAR(1,J)=COVAR(J,J) COVAR(J,J)=0.d0 15 CONTINUE COVAR(LISTA(1),LISTA(1))=SWAP c c...now sort elements into proper order on diagonal c DO 16 J=2,MFIT COVAR(LISTA(J),LISTA(J))=COVAR(1,J) 16 CONTINUE c c...finally fill in above diagonal by symmetry c DO 18 J=2,MA DO 17 I=1,J-1 COVAR(I,J)=COVAR(J,I) 17 CONTINUE 18 CONTINUE c RETURN END C C_____________________________________________________________________________ C SUBROUTINE SUBIF(X,A,Y,DYDA,NA) C C User supplied routine for the Levenberg-Marquardt nonlinear least-squares C fit. For x=X and using the NA constants in A this returns the C ccalulated value of y in Y and partial derivatives of y wrt all the C constants in DYDA. Since the necessary calculations have already been C performed by SETUP only a transfer of data between various arrays is C required here. C C This routine is for subtraction of the IF signal from the interferogram. C The IF (typically 20 MHz) is assumed to be of constant amplitude. C C Y = A(1) * cos ( 2 PI * f_IF * X + A(2) * (2 PI/360)) + A(3) C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=10) parameter (pi=3.141592653589793d0) C REAL*8 X,A(NA),Y,DYDA(NA),aa(na) c common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto c riffr=2.d0*pi*fifmhz phase=a(2)*2.d0*pi/360.d0 Y=a(1)*dcos(riffr*x+phase)+a(3) c c...evaluate the derivatives DYDA numerically c do 1 j=1,na do 2 jj=1,na aa(jj)=a(jj) 2 continue if(j.eq.1)aa(j)=aa(j)+100.d0 if(j.eq.2)aa(j)=aa(j)+2.d0 if(j.eq.3)aa(j)=aa(j)+20.d0 if(j.gt.3)then DYDA(j)=0.d0 goto 1 endif c phase=aa(2)*2.d0*pi/360.d0 YY=aa(1)*dcos(riffr*x+phase)+aa(3) DYDA(j)=(yy-y)/(aa(j)-a(j)) c 1 continue C RETURN END C C_____________________________________________________________________________ C SUBROUTINE SUBLIN(X,A,Y,DYDA,NA) C C User supplied routine for the Levenberg-Marquardt nonlinear least-squares C fit. For x=X and using the NA constants in A this returns the C calculated value of y in Y and partial derivatives of y wrt all the C constants in DYDA. Since the necessary calculations have already been C performed by SETUP only a transfer of data between various arrays is C required here. C C This routine is for subtraction of the selected line from the interferogram, C Initial frequency of the line = A(3) (in MHz) comes from measurement in C the frequency domain. C C Y = A(1) * exp( - X/ A(2) ) * cos [2 PI * A(3) * X + A(4) * (2 PI/360)] + A(5) C IMPLICIT REAL*8 (A-H,O-Z) C parameter (maxpar=10) parameter (pi=3.141592653589793d0) C REAL*8 X,A(NA),Y,DYDA(NA),aa(na) c common /intfre/fifmhz,fifkhz,fifmho,fifkho,ndet,ndeto c riffr=2.d0*pi*a(3) phase=a(4)*2.d0*pi/360.d0 Y=a(1)*dexp(-x/A(2))*dcos(riffr*x+phase)+a(5) c c...evaluate the derivatives DYDA numerically c do 1 j=1,na do 2 jj=1,na aa(jj)=a(jj) 2 continue if(j.eq.1)aa(j)=aa(j)+100.d0 if(j.eq.2)aa(j)=aa(j)+ 0.5d0 c if(j.eq.3)aa(j)=aa(j)+0.001d0 if(j.eq.4)aa(j)=aa(j)+2.d0 c if(j.eq.5)aa(j)=aa(j)+20.d0 if(j.eq.3.or.j.ge.5)then DYDA(j)=0.d0 goto 1 endif c riffr=2.d0*pi*aa(3) phase=aa(4)*2.d0*pi/360.d0 YY=aa(1)*dexp(-x/aa(2))*dcos(riffr*x+phase)+aa(5) DYDA(j)=(yy-y)/(aa(j)-a(j)) c 1 continue C RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________ c c readsp To read a spectral file c summar To display summary information on acquired spectra c marsca To plot and label the frequency marker scale c intcom To compare nine interferograms per screen c startgr To initialise graphics c looksp Main graphics screen c inpout To determine whether .FAR archives are available and c how to input spectra c inpspe To set up a list of files available for input as spectra c FFTEXE To set up the FFT calculation c spectrm Power spectrum estimation c four1 FFT routine c PF Peak measurement c SORTH For sorting input spectra in frequency c baksub Background subtraction c inkey keyboard input c c mrqmin Main routine for Leavenberg-Marquardt fitting method c mrqcof Used by MRMIN to evaluate linearised fitting matrix c gaussj Solution of a set of linear equations using GAUSSJ method c covsrt Deal with covariance matrix from the fit c subif Routine for LM fit defining constant IF wave C_____________________________________________________________________________