c$debug C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C c VIEWF - Display program for individual FTMW spectral files: c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c C Ver 30.I.2001 ---- Zbigniew KISIEL ---- c C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C_________________________/-------------------------------------------------- C C Modification history: C c 16.10.99: Created from FFT5A by chopping out scope and ASCII options, and c adding display features from VIEWM c 20.10.99: Provisional filetering of Doppler doubling trnsferred from F.FOR c 30.01.01: Debugging of O option C_____________________________________________________________________________ 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 The code in the two INCLUDE files clashes with IMPLICIT statements - so C when using MSF5.0 graphics the recomendation is to type everything C explicitly (check with -4Yg compilation switch) C c Uncomment the statements below (and in similar headers in routines c as necessary) MSF5 = MS Fortran 5.0 c PS1 and PS4 = MS Powerstation Fortrans 1.0 and 4.0 c c INCLUDE 'FGRAPH.FI' MSF5+PS1 c INCLUDE 'FLIB.FI' PS1 c USE MSFLIB PS4 C_____________________________________________________________________________ C c---NOTE: MAXSPE should be less than NMAXPT/2, otherwise the EQUIVALENCE c declaration in routine RELTIM may cause problems c PARAMETER (nmaxpt= 8192,maxspe=4000,maxfil=80) c c---------------------------------------------------------------------------- c c idata - compressed and cleaned up interferogram c data - interferogram passed to FFT c p - frequency domain points out of FFT c nposc - number of points in raw interferogram c nrep - repetition count for adding multiple microwave shots c nskips - points to be skipped at the beginning of the interferogram c nskipe - points to be skipped at the end of the interferogram c ncut - number of points in the spectrum out of the FFT c real data(nmaxpt),w1(2*nmaxpt),w2(nmaxpt) real p(nmaxpt) integer idata(maxspe) real*8 fcent,fstep,fmax character filnam*12,fnams(maxfil)*12 common /points/data,npts common /work/w1 common /work1/w2 common /interf/idata,tstep,NREP,NSKIPS,NSKIPE common /spec/p,fstep,ncut,NFFT,NCALL common /sfiles/fnams c C CALL clearscreen($GCLEARSCREEN) WRITE(*,3344) 3344 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | V I E W F - Viewer for FTMW spectral files ', * T79,'|'/ * ' |',76(1H_),'|'/' version 30.I.2001',T64,'Zbigniew KISIEL') c write(*,3345)maxfil,maxspe,nmaxpt 3345 format(1x/' up to',i6,' files can be listed'/ * ' up to',i6,' interferogram points can be read'/ * ' up to',i6,' points for FFT') c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...determine available files c 360 call inpspe(nfil) c c...list available files, columnwise in order of creation c 365 IF(nfil.le.4)then write(*,121)(i,fnams(i),i=1,nfil) goto 123 ENDIF IF(nfil.eq.5)then write(*,121)(i,fnams(i),i=1,4) write(*,124)5,fnams(5) 124 format(61x,i3,' = ',a) goto 123 ENDIF c if((nfil/4)*4.eq.nfil)then j=nfil/4 else j=nfil/4+1 endif jj=j*4-nfil do 120 i=1,j-jj write(*,121)i,fnams(i),i+j,fnams(i+j),i+2*j,fnams(i+2*j), * i+3*j,fnams(i+3*j) 120 continue 121 format(i4,' = ',a,3(2x,i3,' = ',a)) if(j*4.ne.nfil)then nn=j-jj+1 do 122 i=nn,nn+jj-1 write(*,121)i,fnams(i),i+j,fnams(i+j),i+2*j,fnams(i+2*j) 122 continue endif c c...select file c 123 write(*,18) 18 format(1x/' File number (ENTER to exit): ',$) read(*,'(I5)',err=365)nsel if(nsel.lt.0.or.nsel.gt.nfil)goto 365 if(nsel.eq.0)then write(*,'(1x//)') stop endif filnam=fnams(nsel) c call readsp(fcent,vstep,naver,filnam,iread) if(iread.ne.1)goto 365 WRITE(*,'(1X/)') PAUSE "------press ENTER to continue------" c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...FFT c 771 NCALL=0 call FFTEXE c samfre=1./tstep fnyq=samfre/2. fmax=fnyq*1.E-3 c call looksp(fcent,fmax,vstep,filnam) c write(*,367)filnam 367 format(1x/' File just processed: ',a/) goto 365 c stop end c c--------------------------------------------------------------------------- c c SUBROUTINE spctrm(p,m,k,ovrlap) 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= 8192) INTEGER k,m c REAL p(m),w1(2*nmax),w2(nmax) REAL p(nmax),w1(2*nmax),w2(nmax) LOGICAL ovrlap real data(nmax) real*8 fstep common /points/data,npts common /work/w1 common /work1/w2 common /spec/p,fstep,ncut,NFFT,NCALL c INTEGER j,j2,joff,joffn,kk,m4,m43,m44,mm REAL den,facm,facp,sumw,w,window 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(1x/i10,' points used in FFT') c return end c c---------------------------------------------------------------------------- c SUBROUTINE four1(nn,isign) parameter (nmax= 8192) INTEGER isign,nn c REAL data(2*nn) real 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 i,istep,j,m,mmax,n REAL 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 graphicsmode() C C This sets the required graphics mode card. It also C determines the pixel limits on x and y axes (0,maxx), (0,maxy) C INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 C RECORD /rccoord/curpos RECORD /xycoord/ixy RECORD /wxycoord/wxy INTEGER*2 dummy,maxx,maxy,LINOFS,mymode,myrows,mycols integer*4 dummy4,blue,red RECORD /videoconfig/myscreen COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red C C...Determine current graphics mode. C CALL getvideoconfig(myscreen) C C...Set graphics mode C C...this setting was preferred until graphics cards started showing c it as a ribbon narrower than screen height (640*350 and 25*80 chars. EGA) c mymode=$ERESCOLOR c myrows=25 C C...this setting should have widest compatibility (640*480, 30*80 chars. VGA) c (the 16 color mode usually delivers faster scrolling) c mymode=$VRES256COLOR mymode=$VRES16COLOR myrows=30 C C...this is for 800*600, 37*100 chars. SVGA which is available only for MFPS c (the 256 color mode experiences problems with some cards) c mymode=$SRES256COLOR c mymode=$SRES16COLOR c myrows=37 c dummy = setvideomoderows ( mymode, myrows ) myrows=dummy if(myrows.ne.25.and.myrows.ne.30.and.myrows.ne.37)then write(*,'(1x//i5,a,a//)')myrows, * ' <--- this is the number of text lines in this mode', * ' and it is not yet supported' pause endif C IF( dummy .EQ. 0 ) STOP 'Error: cannot set graphics mode' C C...Determine parameters of the graphics mode that is being used C CALL getvideoconfig( myscreen ) maxx = myscreen.numxpixels - 1 maxy = myscreen.numypixels - 1 c linofs=nint(real(maxy)/real(myrows)) mycols = myscreen.numtextcols if(mycols.lt.80)then write(*,'(1x//i5,a,a//)')myrows, * ' <--- this is the number of text columns in this mode', * ' and it is too small' pause endif c c...patch for different use of SETBKCOLOR by MSPS4.0 from MSPS1.0 (and MSF5) c blue=$BLUE red=$RED if($BLUE.ne.#2a0000)then blue=1 red=4 endif c dummy4 = setbkcolor( BLUE ) call clearscreen($GCLEARSCREEN) C RETURN END C C------------------------------------------------------------------------ c subroutine readsp(fcent,vstep,naver,filnam,iread) c PARAMETER (maxspe=4000) real*8 fcent character timdat*20,coment*50,sample*20,INTEXT*30 character filnam*12,errm*40 integer idata(maxspe) common /interf/idata,tstep,NREP,NSKIPS,NSKIPE c OPEN(3,FILE=FILNAM,FORM='BINARY',ERR=503,STATUS='OLD') C errm='reading sample' READ(3,err=503,end=503)sample errm='reading first intext' READ(3,err=503,end=503)intext(1:27) errm='reading nrep,nskips,nskips,fcent' READ(intext,'(3i5,f12.5)',err=503,end=503) * nrep,nskips,nskipe,fcent errm='reading second intext' READ(3,err=503,end=503)intext(1:27) errm='reading tstep,vstep,naver' READ(intext,'(1pe10.3,e10.3,i7)',err=503,end=503) * tstep,vstep,naver errm='reading coment' READ(3,err=503,end=503)coment errm='reading timdat' READ(3,err=503,end=503)timdat c errm='since NREP=0' if(nrep.eq.0)goto 503 errm='since NAVER=0' if(naver.eq.0)goto 503 errm='since FCENT.lt.1000.' if(fcent.lt.1000.)goto 503 errm='since TSTEP.eq.0.0' if(tstep.eq.0.0)goto 503 errm='since VSTEP.eq.0.0' if(vstep.eq.0.0)goto 503 C c...Intensities of spectral data points c errm='reading spectral points: idata()' DO 20 N=1,Nrep if(n.le.maxspe)then READ(3,err=503,end=503)idata(n) else goto 21 endif 20 CONTINUE c 21 if(nrep.gt.maxspe)then nskipe=nskipe-(nrep-maxspe) if(nskipe.lt.1)nskipe=1 nrep=maxspe if(nrep-nskips-nskipe.le.5)then nskipe=5 nskips=5 endif write(*,3470)char(7),maxspe 3470 format(//' ***** Interferogram exceeds program limits and', * ' has been chopped to ', * a,i5,' points') endif C write(*,'(1x//1x,a/)')coment WRITE(*,10)' Sample: ',SAMPLE WRITE(*,10)' Time/Date: ',TIMDAT write(*,11)' No of points: ',nrep write(*,11)' Points skipped at beginning: ',nskips write(*,11)' Points skipped at end: ',nskipe write(*,12)' Microwave frequency (MHz): ',fcent write(*,12)' X-spacing (microseconds): ',tstep/1.E-6 write(*,12)' Y-spacing (Volts) : ',vstep write(*,11)' Number of averages: ',naver c 10 format(1x,2a) 11 format(1x,a,i7) 12 format(1x,a,f20.12) C CLOSE(3) iread=1 return c 503 write(*,504)filnam,errm 504 format(1x/' ***** Cannot open or not spectral file: ',a/ * ' ***** problems ',a/) iread=0 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 PARAMETER (nmaxpt= 8192,maxspe=4000) c real data(nmaxpt),w1(2*nmaxpt),w2(nmaxpt) real p(nmaxpt) integer idata(maxspe) real*8 fstep logical ovrlap common /points/data,npts common /work/w1 common /work1/w2 common /interf/idata,tstep,NREP,NSKIPS,NSKIPE common /spec/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 c 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)npts 56 format(1x/i7,' points readied for FFT') 31 if(NCALL.ne.1)then 33 write(*,30) 30 format(1x/' Specify zero filling parameter n defined as ' * ' 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,char(7) 775 format(1x/' ERROR:',i8,' points needed for FFT but only ', * i6,' dimensioned',a// * ' -----> 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 looksp(fcent,fmax,vstep,filnam) c INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 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 logical*2 true PARAMETER (nmaxpt= 8192,maxspe=4000,maxsmo=199, * nsect=20,true=.true.) c INTEGER*1 MARK [ALLOCATABLE] (:) 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, * hbotsp,htopsp,rinter,cursle,cursri real p(nmaxpt),spol(maxsmo) INTEGER*2 dummy,maxx,maxy,LINOFS,mymode,myrows,mycols,inkey integer*4 dummy4,blue,red character kk,outstr*27,filnam*12 character*80 emplin,lwork1,lwork2,lwork3 integer idata(maxspe),ioldat(maxspe),itemp(maxspe) COMMON /limits/maxx,maxy,LINOFS,curpos,ixy,wxy, * mymode,myrows,mycols,blue,red common /spec/p,fstep,npts,NFFT,NCALL common /interf/idata,tstep,NREP,NSKIPS,NSKIPE common /peak/x(nsect),y(nsect) common /smooth/ioldat,itemp,spol c idefm=0 htcut=0.5 CURINC=1.0 itxt=0 WRITE(emplin,'(80(1H ))') c c...preserve a copy of the interferogram in IOLDAT c do 1099 j=1,nrep ioldat(j)=idata(j) 1099 continue c c. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C C...Start up the graphics C CALL graphicsmode() pixcol=getcolor() dummy4 = setbkcolor( BLUE ) IPIXLN=nint(real(maxy)/real(myrows)) C C...definition of graphics viewport in pixel coordinates C call setviewport(0,0,maxx,maxy) c c...Find Y-limits c 851 fstart=0.d0 fend=500.d0 smin=1.d+20 smax=1.d-20 f=-fstep do 111 i=1,npts f=f+fstep if(f.lt.fstart)goto 111 if(f.gt.fend)goto 112 if(p(i).lt.smin)smin=p(i) if(p(i).gt.smax)smax=p(i) 111 continue 112 smax=1.05*smax smaxol=smax SMIN=0 c top=smax+0.10d0*(smax-smin) bottom=smin-0.1D0*(smax-smin) CDUMP OPEN(7,FILE='DUMP',STATUS='UNKNOWN') CDUMP WRITE(7,'(4F20.8)')SMIN,SMAX,TOP,BOTTOM htopsp=top-(2.*real(ipixln)+2.)/real(maxy)*(top-bottom) itxt=(mycols-80)/2 if(myrows.eq.37)then hbotsp=bottom+(2.*real(ipixln)+1.)/real(maxy)*(top-bottom) else hbotsp=bottom+(1.2*real(ipixln)+1.)/real(maxy)*(top-bottom) endif c c...definition of graphics viewport in floating point coordinates to be used c for plotting C NOTE: for simplicity the whole screen is assigned in this program to the c floating point graphics window. Space corresponding to two top and c one bottom textmode lines is erased by lines carrying descriptive c information. Rectangular boundary frame is drawn by calculating c appropriate floating point coordinates so that vertical pixels just c below and just above descriptive lines are used. c dummy=setwindow(TRUE,fstart,bottom, * fend,top) call clearscreen($GCLEARSCREEN) DUMMY=SETTEXTCOLOR(7) c c...define marker c fmark=(fend-fstart)/2.d0 fincr=(fend-fstart)/maxx if(idefm.eq.0)then $NODEBUG CALL setlinestyle(#9999) c$DEBUG CALL moveto_w(fmark,hbotsp,wxy) dummy=lineto_w(fmark,smax) imsize=imagesize_w(fmark,smax,fmark,hbotsp) ALLOCATE (MARK(imsize)) CALL getimage_w(fmark,smax,fmark,hbotsp,MARK) call putimage_w(fmark,smax,MARK,$GXOR) $NODEBUG CALL setlinestyle(#ffff) c$DEBUG idefm=1 endif c c...Plot c 699 nfirst=0 CALL moveto_w(fstart,hbotsp,wxy) dummy=lineto_w(fend,hbotsp) dummy=lineto_w(fend,htopsp) dummy=lineto_w(fstart,htopsp) dummy=lineto_w(fstart,hbotsp) DO 6 I=1,npts f=(i-1)*fstep if(f.lt.fstart)goto 6 if(f.gt.fend)goto 697 if(nfirst.eq.0)then nfirst=1 RSPEC=p(i) IF(RSPEC.gt.htopsp)rspec=htopsp if(rspec.lt.hbotsp)rspec=hbotsp CALL moveto_w(fstart,dble(rspec),wxy) endif RSPEC=p(i) IF(RSPEC.gt.htopsp)rspec=htopsp if(rspec.lt.hbotsp)rspec=hbotsp dummy=lineto_w(f,RSPEC) 6 CONTINUE c c...marker c 697 call putimage_w(fmark,smax,MARK,$GXOR) c c...information lines c 771 FMIN=FCENT-FMARK/1000.D0 FPLUS=FCENT+FMARK/1000.D0 YVAL=P( NINT(FMARK/FSTEP)+1 ) write(lwork1,'(F7.1,''kHz --> f-'',F11.4,'', f+'',F11.4, * '' MHz'',21X,''Y:'',F8.2)')fmark,FMIN,FPLUS,yval CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork1) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(EMPLIN) CALL settextposition(myrows,int2(itxt+2),curpos) CALL outtext('n='//CHAR(48+NFFT)) write(lwork3,822) smax-smin, fend-fstart 822 format(F12.2,'<- Y_range X_range ->',f8.2,' kHz') CALL settextposition(myrows,int2(itxt+5),curpos) call outtext(lwork3(1:45)) DUMMY=SETTEXTCOLOR(14) CALL settextposition(myrows,int2(itxt+56),curpos) CALL outtext(filnam) DUMMY=SETTEXTCOLOR(7) CALL settextposition(myrows,int2(itxt+72),curpos) CALL outtext('H = help') 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 R - return to initial settings C I - display interferogram c H - help screen c U - ASCII dump of current FFT c O - determine frequency of peak nearest the cursor c 0 - change bisection range for peak measurement c = - central frequency of Doppler doublet (from last two lines c measured with 'O') c P - show the FFT points 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.'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(kk.eq.'U'.or.kk.eq.'u')then open(7,file='f.out',status='unknown') write(7,'(''! Pump frequency ='',f15.6/ * ''! DeltaF power f- f+'')')fcent f=0.d0 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),fmin,fplus f=f+fstep if(f.gt.300000.)goto 33 3 continue 33 close(7) endif c IF(IK.NE.13)GOTO 77 C C...exit C CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(emplin) DUMMY=SETTEXTCOLOR(15) dummy4 = setbkcolor( RED ) WRITE(OUTSTR,'(A)')'ARE YOU SURE ?' CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(outstr(1:14)) CALL settextposition(1,int2(itxt+20),curpos) WRITE(*,'(1X,A1,$)')CHAR(7) 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(7) dummy4 = setbkcolor( BLUE ) CALL settextposition(1,int2(itxt+1),curpos) call outtext(lwork1) GOTO 77 C 915 continue DEALLOCATE(MARK) dummy=setvideomode($DEFAULTMODE) return C C...Shift cursor to the left (with K) C 710 call putimage_w(fmark,smax,MARK,$GXOR) fmark=fmark-8.d0*fincr IF(KK.EQ.'k')fmark=fmark+7.d0*fincr IF(fmark.LT.fstart)fmark=fstart C 719 call putimage_w(fmark,smax,MARK,$GXOR) GOTO 771 C C...Shift cursor to the right (with L) C 711 call putimage_w(fmark,smax,MARK,$GXOR) 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 call putimage_w(fmark,smax,MARK,$GXOR) 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.0.d0)fstart=0.d0 if(fend.gt.fmax)fend=fmax IF(FMARK.LT.fstart)FMARK=Fstart IF(FMARK.GT.fend)FMARK=Fend c 801 fincr=(fend-fstart)/maxx dummy=setwindow(TRUE,fstart,bottom, * fend,top) call clearscreen($GCLEARSCREEN) 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 fstart=0.d0 fend=500.d0 fmark=(fend-fstart)/2.d0 smax=smaxol 802 top= smax+0.10d0*(smax-smin) bottom=SMIN-0.10D0*(smax-smin) CDUMP WRITE(7,'(4F20.8)')SMIN,SMAX,TOP,BOTTOM htopsp=top-(2.*real(ipixln)+2.)/real(maxy)*(top-bottom) if(myrows.eq.37)then hbotsp=bottom+(2.*real(ipixln)+1.)/real(maxy)*(top-bottom) else hbotsp=bottom+(1.2*real(ipixln)+1.)/real(maxy)*(top-bottom) endif goto 801 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 IF(fstart.LT.0.D0)THEN fstart=0.D0 Fchang=fend-(fstart+FRange) CALL settextposition(1,int2(itxt+1),curpos) WRITE(*,101)CHAR(7) 101 FORMAT(1X,A1,$) ENDIF fend=fend-Fchang FMARK=FMARK-Fchang GOTO 801 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 IF(fend.gt.fmax)THEN fend=fmax CALL settextposition(1,int2(itxt+1),curpos) WRITE(*,101)CHAR(7) fstart=fend-frange ENDIF FMARK=FMARK+Fchang GOTO 801 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...Take current marker frequency as measurement of line frequency (with 9) c 1440 xpeak=fmark FLASTM=FNEWM FLASTP=FNEWP flasto=fnewo fNEWM=fcent-dble(xpeak)*0.001d0 FNEWP=FCENT+dble(XPEAK)*0.001d0 fnewo=xpeak write(lwork3,1441)xpeak,'peak',FNEWM,FNEWP 1441 format(F7.1,'kHz <-',a,'-> ',f10.4,' - + ',f10.4) DUMMY=SETTEXTCOLOR(15) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork3) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) DUMMY=SETTEXTCOLOR(7) c DUMMY=SETCOLOR(14) y1=htopsp y2=0.d0 call moveto_w(DBLE(xpeak),DBLE(y1),wxy) dummy=lineto_w(DBLE(xpeak),DBLE(y2)) dummy=setcolor(15) c goto 77 c c...Show the FFT points (with P) C 1300 dely=smax/150. delx=(fend-fstart)/150. DO 1301 i=1,Npts if(p(i).gt.htopsp)goto 1301 f=(i-1)*fstep if(f.lt.fstart)goto 1301 if(f.gt.fend)goto 1302 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...Frequency of peak nearest the cursor (with O) C 1400 FLASTM=FNEWM FLASTP=FNEWP flasto=fnewo call pf(fmark,ypeak,xpeak,errorx,htcut,xbot,fhalm,fhalp) fNEWM=fcent-dble(xpeak)*0.001d0 FNEWP=FCENT+dble(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 write(lwork3,1410)xpeak,'peak',FNEWM,FNEWP,errorx*1000.,ypeak 1410 format(F7.1,'kHz <-',a,'-> ',f10.4,' - + ',f10.4, * ' (Er.fit=',f6.1,') Y:',f8.2) DUMMY=SETTEXTCOLOR(15) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork3) WRITE(lwork3,'(67x,''FWHH:'',F8.2)')fhalp-fhalm CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(lwork3) DUMMY=SETTEXTCOLOR(7) c c...bisection points dely=smax/150. delx=(fend-fstart)/150. DUMMY=SETCOLOR(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(12) y2=0.4*ypeak if(y2.lt.top)then call moveto_w(DBLE(xpeak),DBLE(ypeak),wxy) dummy=lineto_w(DBLE(xbot),DBLE(y2)) endif c c...vertical bar DUMMY=SETCOLOR(14) y1=ypeak y2=0.d0 if(y1.lt.top)then call moveto_w(DBLE(xpeak),DBLE(y1),wxy) dummy=lineto_w(DBLE(xpeak),DBLE(y2)) endif c c...horizontal bars at 0, 0.5 and 1.0 Ymax x1=xpeak-5.*delx x2=xpeak+5.*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(x2,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(15) c goto 77 c c...Mean frequency of last two measured peaks (for Doppler doublets), with = C 1430 fmin=(flastm+fnewm)*0.5d0 fplus=(flastp+fnewp)*0.5d0 foffs=(flasto+fnewo)*0.5d0 write(lwork3,1411)foffs,'mean',fmin,fplus,abs(flasto-fnewo) 1411 format(F7.1,'kHz <-',a,'-> ',f10.4,' - + ',f10.4,' splitting', * f7.2,' kHz') DUMMY=SETTEXTCOLOR(11) CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(lwork3) DUMMY=SETTEXTCOLOR(7) CALL settextposition(2,int2(itxt+1),curpos) CALL outtext(emplin) 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, * fend,top) call clearscreen($GVIEWPORT) GOTO 699 c c . . . . . . . . . . . . . . . . . . . . . . . . c c...display the interferogram (with I) and allow selection of discarded c points with scrollable cursors c 830 ICHANG=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=setwindow(TRUE,1.0d0,dble(botlim), * DBLE(nrep),dble(toplim)) call clearscreen($GVIEWPORT) 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) dummy=setcolor(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 dummy=setcolor(15) DO 86 I=nskips+1,nrep-nskipe Rinter=DBLE(Idata(I)) dummy=lineto_w(DBLE(I),Rinter) 86 CONTINUE dummy=setcolor(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 dummy=setcolor(15) c CURSLE=DBLE(NSKIPS) CURSRI=DBLE(NREP-NSKIPE) call putimage_w(CURSLE,dble(maxi),MARK,$GXOR) call putimage_w(CURSRI,dble(maxi),MARK,$GXOR) C WRITE(LWORK2,'(2A)') * 'A,S <-cursors-> K,L + - N B R O G ', * 'U T H=help ENTER to exit' CALL settextposition(1,int2(itxt+1),curpos) call outtext(lwork2) write(lwork3,82) (maxi-mini)*vstep*1000,tstep*1.d6 82 format(' Yrange ',f9.2,'mV ',F5.2,'us/pt') CALL settextposition(myrows,int2(itxt+1),curpos) call outtext(lwork3(1:38)) CALL settextposition(myrows,int2(itxt+52),curpos) WRITE(OUTSTR,1550)nint(cursle),nint(cursri), * (nint(cursri)-nint(cursle))*tstep*1.e6 1550 format('in use:',i5,',',i4,' ->',F5.1,'us') call outtext(outstr) 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 U c - relaxation time T 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.'T'.OR.KK.EQ.'t')GOTO 1600 if(kk.eq.'G'.or.kk.eq.'g')goto 1700 if(kk.eq.'O'.or.kk.eq.'o')goto 1800 if(kk.eq.'H'.or.kk.eq.'h')goto 1900 IF(KK.EQ.'+'.OR.KK.EQ.'=')CURINC=CURINC*5.0 IF(KK.EQ.'-'.OR.KK.EQ.'_')CURINC=CURINC*0.2 if(kk.eq.'U'.or.kk.eq.'u')then open(7,file='t.out',status='unknown') write(7,'(''! Pump frequency/MHz ='',f15.6)')fcent write(7,'(''! Point spacing/microsec ='',f15.6/''!''/ * ''! time voltage/mV)'')')tstep/1.e-6 do 14 j=1,nrep write(7,15)real(j)*tstep/1.E-6,real(idata(j))*1000.d0*vstep 14 continue 15 format(f5.1,f10.3) close(7) 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(myrows,int2(itxt+2),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 835 ENDIF c C...Shift RIGHT cursor to the left (with K) C 870 call putimage_w(CURSRI,dble(maxi),MARK,$GXOR) CURSRI=CURSRI-8.d0*CURINC IF(KK.EQ.'k')CURSRI=CURSRI+7.d0*CURINC IF(CURSRI.LE.CURSLE)CURSRI=CURSLE+CURINC C 872 call putimage_w(CURSRI,dble(maxi),MARK,$GXOR) ICHANG=1 CALL settextposition(myrows,int2(itxt+52),curpos) WRITE(OUTSTR,1550)nint(cursle),nint(cursri), * (nint(cursri)-nint(cursle))*tstep*1.e6 call outtext(outstr) GOTO 834 C C...Shift RIGHT cursor to the right (with L) C 875 call putimage_w(CURSRI,dble(maxi),MARK,$GXOR) CURSRI=CURSRI+8.d0*CURINC IF(KK.EQ.'l')CURSRI=CURSRI-7.d0*CURINC IF(CURSRI.GT.dble(NREP))CURSRI=dble(NREP) GOTO 872 c C...Shift LEFT cursor to the left (with A) C 880 call putimage_w(CURSLE,dble(maxi),MARK,$GXOR) CURSLE=CURSLE-8.d0*CURINC IF(KK.EQ.'a')CURSLE=CURSLE+7.d0*CURINC IF(CURSLE.LT.0.d0)CURSLE=0.d0 C 882 call putimage_w(CURSLE,dble(maxi),MARK,$GXOR) ICHANG=1 CALL settextposition(myrows,int2(itxt+52),curpos) WRITE(OUTSTR,1550)nint(cursle),nint(cursri), * (nint(cursri)-nint(cursle))*tstep*1.e6 call outtext(outstr) GOTO 834 c C...Shift LEFT cursor to the right (with S) C 885 call putimage_w(cursle,dble(maxi),MARK,$GXOR) cursle=cursle+8.d0*CURINC IF(KK.EQ.'s')cursle=cursle-7.d0*CURINC IF(cursle.ge.cursri)cursle=cursri-1.d0*CURINC GOTO 882 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//12x,'Number of points in the interferogram: ',i5/12x, * 'Number of points (>3 and odd) in smoothing interval .... ',$) READ(*,'(i5)',ERR=1101)NSPT IF(NSPT.LE.3.OR.NSPT.GT.maxsmo)GOTO 1101 IF((NSPT/2)*2.EQ.NSPT)GOTO 1101 WRITE(*,'(1X//'' S U B T R A C T I N G''//)') C call baksub(nspt,1) c ICHANG=1 GOTO 831 C C...Compensate for rotational relaxation (with N) C 1500 CALL clearscreen($GCLEARSCREEN) 1501 WRITE(*,1502)nrep 1502 FORMAT(1X//12x,'Number of points in the interferogram: ',i5/ * 12x,'Number of points for halfdecay .... ',$) READ(*,'(i5)',ERR=1501)NHALF IF(NHALF.LE.10.AND.NHALF.GE.-10)GOTO 1501 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...Relaxation time (with T) C 1600 do 1601 j=1,nrep idata(j)=iabs(idata(j)) 1601 continue c 1603 CALL settextposition(1,int2(itxt+1),curpos) CALL outtext(emplin(1:79)) CALL settextposition(myrows,int2(itxt+1),curpos) CALL outtext(emplin(1:79)) CALL settextposition(1,int2(itxt+1),curpos) write(*,'(1H+, * '' No of points (>3 and odd) for smoothing: '',$)') READ(*,'(i5)',ERR=1603)NSPT IF(NSPT.LE.3.OR.NSPT.GT.maxsmo)GOTO 1603 IF((NSPT/2)*2.EQ.NSPT)GOTO 1603 CALL settextposition(myrows/2+5,int2(itxt+1),curpos) c call baksub(nspt,0) c c...fit y=a exp(-bt), time for halfdecay=-exp(0.5)/b call reltim(a,b,erb) c dummy=setcolor(int2(14)) rinter=idata(nskips) if(idata(nskips).gt.maxi)rinter=maxi if(idata(nskips).lt.mini)rinter=mini CALL moveto_w(DBLE(nskips),rinter,wxy) DO 1606 I=nskips+1,nrep-nskipe Rinter=DBLE(Idata(I)) dummy=lineto_w(DBLE(I),Rinter) 1606 CONTINUE c dummy=setcolor(int2(11)) t=tstep/1.E-6 rinter=a*exp(-b*t) if(idata(nskips).gt.maxi)rinter=maxi if(idata(nskips).lt.mini)rinter=mini CALL moveto_w(DBLE(nskips),rinter,wxy) DO 1608 I=nskips+1,nrep-nskipe t=(i-nskips+1)*tstep/1.E-6 rinter=a*exp(-b*t) dummy=lineto_w(DBLE(I),Rinter) 1608 CONTINUE c thalf=0.69314718/b ethalf=0.69314718*erb/(b*b) write(*,1605)a,b,thalf,ethalf 1605 format(' A =',f15.7/ * ' B =',f15.7,5x, * ' t(1/2) =',f9.5,' +-',f8.5,' microseconds') CALL settextposition(myrows,int2(itxt+1),curpos) CALL outtext(' Press E N T E R to continue ') c 1607 IK=INKEY(N) if(ik.ne.13)goto 1607 c ICHANG=1 GOTO 1200 c c...filter-out Doppler doubling (with G) c 1700 call decpar(fdoubl,phase,npatch,nsmopt, * maxi,mini,ichang) goto 831 c c...Replace inteferogram in IOLDAT with the current interferogram (with O) c 1800 do 1801 j=1,nrep ioldat(j)=idata(j) 1801 continue ichang=0 goto 831 c c...Restore original interferogram (with R) c 1200 do 1201 j=1,nrep idata(j)=ioldat(j) 1201 continue ichang=1 goto 831 c c...The interferogram help screen c 1900 CALL clearscreen($GCLEARSCREEN) write(*,'(1x)') WRITE(*,1901) 1901 FORMAT( ' COMMANDS ACTIVE IN THE INTERFEROGRAM DISPLAY:'/ * ' ---------------------------------------------'// * ' A/S - shift the left hand cursor'/ * ' K/L - shift the right hand cursor'/ * ' caps on/off - fast/slow change in the above'/ * ' + - - change the increment of cursor shift'// * ' B - subtract the low frequency baseline'/ * ' R - restore the original interferogram'/ * ' O - replace original interferogram by current one'/ * ' U - ASCII dump of current interferogram'// * ' T - determine relaxation time'/ * ' G - filter out Doppler doubling'// * ' Exit back to the FFT screen by pressing ENTER'//) 1902 IK=INKEY(J) IF(IK.NE.13)GOTO 1902 dummy=setwindow(TRUE,fstart,bottom, * fend,top) ichang=1 GOTO 831 C c . . . . . . . . . . . . . . . . . . . . . . . . c c...display the main help screen (with H) c 840 CONTINUE CALL clearscreen($GCLEARSCREEN) write(*,'(1x)') WRITE(*,232) 232 FORMAT( ' SUMMARY OF COMMANDS ACTIVE IN FFT MODE:'/ * ' ---------------------------------------'// * ' 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 - show the interferogram'/ * ' N - change FFT zero filling parameter n'/ * ' P - show the FFT points'/ * ' R - rescale spectrum to initial conditions'/ * ' U - ASCII dump of current FFT'/ * ' 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'// * ' Exit from the graphics is by ENTER followed by Y'/) WRITE(*,106) 106 FORMAT(40X, * 'ENTER returns to the spectral display ',$) 107 IK=INKEY(J) IF(IK.NE.13)GOTO 107 835 dummy=setwindow(TRUE,fstart,bottom, * fend,top) call clearscreen($GVIEWPORT) GOTO 699 c c...change FFT zero filling parameter n c 850 CALL clearscreen($GCLEARSCREEN) call FFTEXE GOTO 851 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= 8192,nsect=20) c real p(nmaxpt),fr(nsect,2) real*8 fstep,fmark REAL*8 SUMX,SUMY,SUMXY,SUMX2,SUMY2,CXX,CXY,CYY,RN,A0,A1, * xoffs,yoffs c common /spec/p,fstep,npts,NFFT,NCALL common /peak/x(nsect),y(nsect) 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 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 15 continue c RETURN END C_____________________________________________________________________________ c cc integer*2 function INKEY(N2) msps 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 aniother call to GETCHARQQ is required to c get the extended code of the character c cc INCLUDE 'FLIB.FD' PS1 c USE MSFLIB PS4 c cc INTEGER*2 IK msps cc CHARACTER*1 KK msps c cc KK=GETCHARQQ() msps cc IK=ICHAR(KK) msps cc IF(IK.EQ.0 .OR. IK.EQ.224 ) THEN msps cc KK=GETCHARQQ() msps cc IK=-ICHAR(KK) msps cc ENDIF msps cc n2=ik msps cc INKEY=IK msps cc END msps C C------------------------------------------------------------------------ interface to integer*2 function system [C] (string[reference]) msf5 character*1 string msf5 end msf5 C------------------------------------------------------------------------ c subroutine inpspe(nfil) c c Routine to set up a list of files available for input as spectra c c c INCLUDE 'FLIB.FD' PS1 c USE MSFLIB PS4 c integer*2 system msf5 logical*4 fsys parameter (maxfil=80) character line*20,fnams(maxfil)*12 common /sfiles/fnams c c c...save the current directory to file c write(*,'(1x/)') c fsys=systemqq('dir /Od>spec.lst') msps ifsys= system('dir /Od>spec.lst'c) msf5 fsys=.true. msf5 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 nfil=0 7 read(2,'(a)',end=9)line 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')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(1:6).ne.'INTERF'.and. * (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.maxfil)then write(*,'(1x//'' ***** The current limit of'',i3,a * '' on the number of files has been reached''/ * '' ***** NO MORE WILL BE READ IN''/)') * maxfil,char(7) goto 9 endif goto 7 9 close(2) c fsys=systemqq('del spec.lst') msps ifsys= system('del spec.lst'C) msf5 c c...options c 19 if(nfil.eq.0)then write(*,100) 100 format(1x//' S O R R Y:'//5x, * ' This directory appears to contain NO FTMW spectral files'//) stop endif c return end c C_____________________________________________________________________________ c subroutine baksub(nspt,nsub) c c NSPT - number of points in the smoothing interval c NSUB = 1 subtract the smoothed spectrum from the original c = 0 do not subtract c PARAMETER (nmaxpt= 8192,maxspe=4000,maxsmo=199) c integer idata(maxspe),ioldat(maxspe),itemp(maxspe) real spol(maxsmo) common /interf/idata,tstep,nrep,nskips,nskipe 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 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=-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 the smoothed background in IDATA from the spectrum in IOLDAT c if(nsub.eq.1)then do 1106 j=1,nrep idata(j)=ioldat(j)-idata(j) 1106 continue endif c return end c C_____________________________________________________________________________ c subroutine reltim(a,b,erb) c c The expression c y = A exp(-B t) c is fitted to the data in IDATA via least squares straight line fit of c ln(y) = ln(A) - B t c c---NOTE: MAXSPE should be less than NMAXPT/2, otherwise the EQUIVALENCE c declaration may cause problems c PARAMETER (nmaxpt= 8192,maxspe=4000,maxsmo=199) c integer idata(maxspe) common /work/w1(2*nmaxpt) real*8 x(maxspe),y(maxspe) REAL*8 SUMX,SUMY,SUMXY,SUMX2,SUMY2,CXX,CXY,CYY,RN,A0,A1 equivalence (w1(1),x(1)),(w1(nmaxpt+1),y(1)) common /interf/idata,tstep,nrep,nskips,nskipe c c...set up X and Y (TSTEP is in seconds, while X is set up in microsec) c npts=0 do 1 i=1,nrep if(i.lt.nskips)goto 1 if(i.gt.(nrep-nskipe))goto 2 npts=npts+1 x(npts)=npts*tstep/1.E-6 if(idata(i).gt.0)then y(npts)=dlog(dble(idata(i))) else y(npts)=dlog(dble(idata(1))) endif 1 continue 2 continue c open(2,file='xy.dat',status='unknown') c do 300 i=1,npts c write(2,*)x(i),y(i),idata(i+nskips-1) c300 continue c close(2) c c...least squares c SUMX=0.D0 SUMY=0.D0 SUMXY=0.D0 SUMX2=0.D0 SUMY2=0.D0 DO 7 I=1,npts SUMy=SUMy+Y(I) SUMx=SUMx+X(I) SUMXY=SUMXY+X(I)*Y(I) SUMy2=SUMy2+Y(I)**2 SUMx2=SUMx2+X(I)**2 7 CONTINUE C C...coefficients RN=NPTS CXX=SUMX2-SUMX*SUMX/RN CXY=SUMXY-SUMX*SUMY/RN A1=CXY/CXX A0=(SUMY-A1*SUMX)/RN C C...error CYY=SUMY2-SUMY*SUMY/RN ERA1S=((CYY/CXX)-(CXY/CXX)**2)/(RN-2.D0) ERA0S=SUMX2*ERA1S/RN c b=-a1 ERb=dsqrt(dble(era1s)) a=exp(a0) c write(*,10)npts 10 format(' Results of fitting y = A exp (-Bt) to',i4,' points:'/) c return end c C------------------------------------------------------------------------ c subroutine decpar(fdoubl,phase,npatch,nsmopt, * maxi,mini,ichang) c INCLUDE 'FGRAPH.FD' MSF5+PS1 c USE MSFLIB PS4 C RECORD /wxycoord/wxy c PARAMETER (Nmaxpt= 8192,maxspe=4000,maxsmo=199) c real spol(maxsmo) integer*2 dummy,inkey character*1 kk real*8 rinter integer idata(maxspe),ioldat(maxspe),itemp(maxspe) common /interf/idata,tstep,NREP,NSKIPS,NSKIPE common /smooth/ioldat,itemp,spol c CALL clearscreen($GCLEARSCREEN) write(*,'(a)')' ' if(fdoubl.eq.0.0)fdoubl=nrep/2 1601 WRITE(*,'('' NOTE: doubling period = 1/(Doppler doubling *'', * '' point spacing)''//)') WRITE(*,1602)' Doubling period /points',fdoubl 1602 FORMAT(1x,a,' = ',f10.2,' New = ',$) read(*,1603,err=1601)rtemp 1603 format(f20.10) if(rtemp.ne.0.)then if(rtemp.gt.20.0.and.rtemp.lt.real(nrep))then fdoubl=rtemp else goto 1601 endif endif c 1605 WRITE(*,1602)' Node position (phase) /points',phase read(*,1603,err=1605)rtemp if(rtemp.ne.0.)then if(rtemp.gt.0.0.and.rtemp.lt.real(nrep))then phase=rtemp else goto 1605 endif endif c 1620 if(npatch.eq.0)npatch=5 WRITE(*,1621)' Patch length /odd points',npatch 1621 FORMAT(1x,a' = ',i7,' New = ',$) read(*,'(i10)',err=1620)ntemp if(ntemp.ne.0)then if(ntemp.ge.1.and.(ntemp-(ntemp/2)*2).eq.1)then npatch=ntemp else goto 1620 endif endif c 1630 if(nsmopt.eq.0)nsmopt=5 WRITE(*,1621)' Smoothing window /odd points',nsmopt read(*,'(i10)',err=1630)ntemp if(ntemp.ne.0)then if(ntemp.ge.3.and.(ntemp-(ntemp/2)*2).eq.1)then if(ntemp.gt.49)goto 1630 nsmopt=ntemp else goto 1630 endif endif c c...calculate the filtering sine wave c freqsn=3.14159/fdoubl do 1607 j=1,nrep itemp(j)=sin(freqsn*j-freqsn*phase)*1000000. 1607 continue c c...plot filtering sine wave c dummy=setcolor(14) rinter=real(itemp(1))/1000000.*0.25*(maxi-mini) if(rinter.gt.maxi)rinter=maxi if(rinter.lt.mini)rinter=mini CALL moveto_w(DBLE(1),rinter,wxy) DO 1608 I=1,nrep Rinter=DBLE(Itemp(I))/1000000.*0.25*(maxi-mini) if(rinter.gt.maxi)rinter=maxi if(rinter.lt.mini)rinter=mini dummy=lineto_w(DBLE(I),Rinter) 1608 CONTINUE call moveto_w(dble(1),0.0d0,wxy) dummy=lineto_w(dble(nrep),0.0d0) c c...plot data c dummy=setcolor(15) rinter=idata(1) if(idata(1).gt.maxi)rinter=maxi if(idata(1).lt.mini)rinter=mini CALL moveto_w(DBLE(1),rinter,wxy) DO 1609 I=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) 1609 CONTINUE c c...calculate and plot filtered interferogram: c dummy=setcolor(12) rinter=real(itemp(1))/1000000. if(abs(rinter).lt.0.2)rinter=sign(0.2,rinter) rinter=idata(1)/rinter itemp(1)=rinter if(rinter.gt.real(maxi))rinter=maxi if(rinter.lt.real(mini))rinter=mini CALL moveto_w(DBLE(1),rinter,wxy) DO 1611 I=1,nrep rinter=real(itemp(i))/1000000. if(abs(rinter).lt.0.2)rinter=sign(0.2,rinter) itemp(i)=idata(i)/rinter rinter=itemp(i) if(rinter.gt.real(maxi))rinter=maxi if(rinter.lt.real(mini))rinter=mini dummy=lineto_w(DBLE(I),Rinter) 1611 CONTINUE c c...show patches c dummy=setcolor(11) rtemp=phase-3.0*fdoubl 1614 rtemp=rtemp+fdoubl i=nint(rtemp) ntemp=npatch/2+1 if(i.lt.ntemp+1)goto 1614 if(i.gt.nrep-ntemp)goto 1610 call moveto_w(dble(i-ntemp),dble(itemp(i-ntemp)),wxy) ik=0 do 1616 j=i-ntemp+1,i+ntemp-1 ik=ik+1 rinter=itemp(i-ntemp)+ * real(itemp(i+ntemp)-itemp(i-ntemp))*ik/(2.0*real(ntemp)) dummy=lineto_w(dble(j),rinter) 1616 continue dummy=lineto_w(dble(i+ntemp),dble(itemp(i+ntemp))) goto 1614 c c c...reject filtering with ENTER or accept with 'A' and transfer to IDATA c 1610 ik=inkey(n) kk=char(ik) if(kk.eq.'A'.or.kk.eq.'a')then do 1612 i=1,nrep idata(i)=itemp(i) 1612 continue rtemp=phase-3.0*fdoubl 1625 rtemp=rtemp+fdoubl i=nint(rtemp) if(i.lt.ntemp+1)goto 1625 if(i.gt.nrep-ntemp)goto 1615 ik=0 do 1618 j=i-ntemp+1,i+ntemp-1 ik=ik+1 idata(j)=idata(i-ntemp)+ * real(idata(i+ntemp)-idata(i-ntemp))*ik/(2.0*real(ntemp)) 1618 continue goto 1625 c 1615 ichang=1 c c...smooth jagged edges of patches c c...set up coefficients in smoothing polynomial M=NSmoPT/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,NSmoPT IS=j-M-1 SPOL(j)=T1*(T2-5*IS*IS) 1103 CONTINUE C c...Smooth three times 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,NSmoPT IS=J-M-1 II=I+IS IF(II.LT.1)II=-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 return endif c if(ik.ne.13)goto 1610 ichang=0 c return end c C------------------------------------------------------------------------ C_____________________________________________________________________________ C_____________________________________________________________________________