C----------------------------------------------------------------------------- c c F R E C A L - Frequency correction in broadband spectra C c Carried out by means of a spline function fitted C to previously measured calibration points. C C C version 8.VI.2009 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C C 8.06.09: Creation from IC C C----------------------------------------------------------------------------- C C TYPICAL OPERATION: C C 1/ read the calibration data in .FRE standard determined in AABS by C using the F8 option in ASCP_L C 2/ determine spline function describing the data and produce output C for checking its quality C 3/ read the spectrum C 4/ recalibrate spectrum using the spline function C 5/ save spectrum C C ALTERNATIVELY: C C If a separate reference spectrum has been recorded at C the same nominal frequency points then a spline function determined C from that spectrum can be read in and used C C Much intermediate output is produced C C----------------------------------------------------------------------------- C C>>> INPUT files, where NAME is the generic name specified inn the input C dialogue: C C 1 / name.spe = spectrum to be rescaled in IFPAN binary standard C 2/ FRECAL.INP = control file for customising interpolation and C smoothing in determining the calibration function C C and either: C C 3a/ name.fre = file of measurements of frequencies and intensities of c lines in the spectrum to be rescaled with added c calculated intensities c (in extended FRE standard of SVIEW peakfinder, entries c generated using the F8 option of SVIEW_L) c or: c C 3b/ name_spline.fnc C C C>>> OUTPUT files: C C 1/ name_calpts.out = x,y listing of frequencies and scaling coefficients i.e. C values of obs/calc intensity scaled around the mean value C (set to unity) c C 2/ name_spline.out = x,y listing of the calculated spline function at c 10 times the point spacing in the spectrum c c 3/ name_spline.fnc = parameters defining the spline function describing c the data c C 4/ name_frecal.spe = 2 column ASCII output of the frequency calibrated C spectrum. This is output for all the points for which C Y-values were defined and can be made equidistant C in frequency with SVIEW_L C C C///////////////////////////////////////////////////////////////////////////// C MODULE globdat integer*4 maxspl PARAMETER (MAXSPL=50000) c real*8 X(maxspl),Y(maxspl),A(maxspl),B(maxspl),C(maxspl), * D(maxspl),GRA,CA,GRB,CB real*8 scalc(maxspl),xold(maxspl),yold(maxspl) integer*4 nspl END MODULE globdat C C///////////////////////////////////////////////////////////////////////////// C MODULE specdat integer*4 maxpts parameter (maxpts=5000000) integer*4 ispec(maxpts) real*8 FREQS(MAXPTS),freqsc(MAXPTS),ftemp(maxpts) REAL*8 FSTART,FEND,FINCR,baslev INTEGER*4 NPTS,ISMALL,ILARGE REAL*4 VKMST,VKMEND,GRID,SAMPRE,GAIN,TIMEC,PHASE, 1 PPS,FRMOD,FRAMPL INTEGER*2 IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC CHARACTER*72 COMENT CHARACTER*6 LAMP,SCANSP CHARACTER*20 SAMPLE CHARACTER*70 FILNAM END MODULE specdat C C///////////////////////////////////////////////////////////////////////////// C C use globdat use specdat C implicit real*8 (a-h,o-z) parameter (maxsmo=10001) character*100 line,filgle character*100 fgen,finsp,ffre,fcalpts,fsplpts,foutsp,fsplfnc real*8 spol(maxsmo) c c WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | FRECAL - Frequency Correction of Broadband Spectra ', * T79,'|'/ * ' |',76(1H_),'|'/' version 8.VI.2009',T64,'Zbigniew KISIEL'/) c c...generic name for files c 100 write(*,101) 101 format(1x//' Generic name to be used for input/output files: ',$) read(*,'(a)')fgen c c...read the control file FRECAL.INP with tuning parameters c c NINTP - the number of linearly interpolated points between data c points defining the spline function (a significant number c may be required for a strongly oscillating calibration) c FIMULT - the multiplier of the point spacing in the spectrum c defining the point spacing to be used for generating the c spline function trace written to NAME_SPLINE.OUT c fimult=10.d0 c open(1,file='frecal.inp',err=200,status='OLD') 27 read(1,'(a)',err=200)line if(line(1:1).eq.'!')goto 27 read(line(44:),*,err=200)nintp if(nintp.lt.0.or.nintp.gt.50)nintp=0 c read(1,'(a)')line read(line(44:),*,err=200)fsmooth if(fsmooth.lt.0.d0.or.fsmooth.gt.50000.d0)fsmooth=0.d0 close(1) goto 210 200 write(*,'(1x//'' ***** ERROR: problem reading FRECAL>INP'',a//)') stop c c...Files for data input/output C 210 finsp=fgen(1:len_trim(fgen))//'.spe' ffre=fgen(1:len_trim(fgen))//'.fre' fcalpts=fgen(1:len_trim(fgen))//'_calpts.out' fsplpts=fgen(1:len_trim(fgen))//'_spline.out' foutsp=fgen(1:len_trim(fgen))//'_frecal.spe' fsplfnc=fgen(1:len_trim(fgen))//'_spline.fnc' filgle=fgen(1:len_trim(fgen))//'.gle' c C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...read the experimental spectrum C call readsp(finsp) <----- C C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C 52 write(*,51) 51 format(//' Specify the source of the calibration data:'// * 19x,'0 = points'/ * 19x,'1 = spline function'// * 14x,'.... ',$) read(*,'(i5)',err=52)inflag c c...Input a previously evaluated spline function c if(inflag.eq.1)then OPEN(3,FILE=fsplfnc,STATUS='unknown',FORM='BINARY',ERR=55) line='Spline function from points in '//ffre(1:len_trim(ffre)) read(3)line(1:100) read(3)nspl do 50 n=1,nspl read(3)x(n),a(n),b(n),c(n),d(n) 50 continue read(3)gra,ca,grb,cb read(3)splst,splinc if(splst.ne.fstart.or.splinc.ne.fincr)then write(*,54)fstart,fincr,splst,splinc 54 format(1x/' ERROR: correct parameter mismatch'// * 10x,'spectrum start:',F15.6/ * 10x,'spectrum step: ',F15.6// * 10x,' spline start:',F15.6/ * 10x,' spline step: ',F15.6//) close(3) stop endif close(3) goto 53 55 write(*,56)fsplfnc(1:len_trim(fsplfnc)) 56 format(1x//' ERROR: cannot read or open ',a//) stop endif c C C...Input the calibration data, which have to be in the version of the F8 standard C of ASCP_L from 2008: c - the first four values are measured from the spectrum in SVIEW_L and c are in FRE standard c - the remaining values are from ASCP_L and are for prediction corresponding c to the measured line c - last number is freq(obs-calc) c c 182705.89149 +- 0.02849 26669 6.77413 182706.04760 3.0648E-04 1.00000 A 14 1 13 <- 13 2 12 -0.15611 c C write(*,'(1x/'' Reading: '',A/)')ffre(1:len_trim(ffre)) open(3,file=ffre,status='old',err=12) C nlin=0 nblend=1 flast=0.d0 c c...Input lines c c Y-axis points are calculated-observed frequency differences c 1 read(3,'(a)',end=7)line c read(line,2,err=1)fobs,erfobs,intobs,fwhh,fcalc 2 format(f14.5,3x,f8.5,i8,f11.5,f16.5) if(fobs.eq.0.d0)goto 1 if(intobs.eq.0)goto 1 if(fcalc.eq.0.d0)goto 1 if(fwhh.eq.0.0d0)goto 1 c c...deal with blends (so far treated as a simple frequency mean) c if(fobs.eq.flast)then if(nlin.eq.0)goto 1 nblend=nblend+1 y(nlin)=(y(nlin)*(nblend-1)+(fcalc-fobs))/real(nblend) write(*,'('' Blend at:'',F15.5)')x(nlin) else nlin=nlin+1 x(nlin)=fobs y(nlin)=fcalc-fobs flast=fobs nblend=1 endif c goto 1 c 7 close(3) nspl=nlin write(*,'(1x/i9,'' calibration points have been read'')')nspl c c...Convert frequency to c scaled frequency X/FSTART to improve conditioning c do 40 n=1,nspl x(n)=x(n)/fstart 40 continue C C...Sort data points in order of increasing X C DO 201 I=1,NSPL-1 J=I 206 J=J+1 IF(X(J)-X(I).lt.0.d0)THEN xx=X(I) X(I)=X(J) X(J)=xx yy=Y(I) Y(I)=Y(J) Y(J)=yy endif IF(J.EQ.NSPL)GOTO 201 GOTO 206 201 CONTINUE c c...write out scaled input points and convert frequency to c scaled frequency X/FSTART to improve conditioning c open(2,file=fcalpts,status='unknown',err=12) write(2,59)ffre(1:len_trim(ffre)) 59 format('!'/'! Sorted calibration points from: ',a/'!') do 58 n=1,nspl xx=x(n)*fstart write(2,22)xx,y(n) if(n.gt.1)then if((x(n)-x(n-1))*fstart.lt.10.d0*fincr)write(*,60)xx endif 60 format(' WARNING: calibration points may be too close at',F15.5) 58 continue close(2) c c...Interpolate NINTP points between existing points in order to reduce c excursions of spline function and then insert the interpolated points c fraint=1.d0/real(nintp+1) c do 63 n=1,nspl xold(n)=x(n) yold(n)=y(n) 63 continue c nn=1 do 61 n=2,nspl do 62 nnn=1,nintp nn=nn+1 x(nn)=xold(n-1)+nnn*fraint*(xold(n)-xold(n-1)) y(nn)=yold(n-1)+nnn*fraint*(yold(n)-yold(n-1)) 62 continue nn=nn+1 x(nn)=xold(n) y(nn)=yold(n) if(nn.gt.maxspl-nintp)then write(*,'(1x//'' ERROR: too few spline points - increase'', * '' MAXSPL''//)') stop endif 61 continue nspl=nn write(*,'(i9,'' points readied for spline calculation''/)')nspl C C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c C...evaluate spline function C call spline(fstart) <----- c C...create binary file of function C OPEN(2,FILE=fsplfnc,STATUS='unknown',FORM='BINARY',ERR=12) line='Spline function from points in '//ffre(1:len_trim(ffre)) write(2)line(1:100) write(2)nspl do 57 n=1,nspl write(2)x(n),a(n),b(n),c(n),d(n) 57 continue write(2)gra,ca,grb,cb write(2)fstart,fincr close(2) C C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...Apply spline frequency correction by evaluating new c frequency at each spectral point c 53 write(*,31) 31 format(1x/' Determining frequency correction to point:'/1x,$) c do 105 n=1,npts call valy(yy,freqs(n)/fstart,nfunct) <----- freqsc(n)=yy if((n/100000)*100000.eq.n)write(*,'(1H+,i8,$)')n 105 continue write(*,'(1H+,i8/)')n c c...Smooth the frequency correction if so specified c if(fsmooth.eq.0.d0)goto 106 nspt=fsmooth/fincr if((nspt/2)*2.eq.nspt)nspt=nspt+1 if(nspt.gt.maxsmo)nspt=maxsmo write(*,102)nspt 102 format(1x/' Smoothing with',i5,' point polynomial'/) 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 C 1200 M=NSPT/2 MM=NSPT 1201 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,MM IS=j-M-1 SPOL(j)=T1*(T2-5*IS*IS) 1103 CONTINUE c c...straightforward smoothing c do 2105 j=1,npts ftemp(j)=freqsc(j) 2105 continue c 2545 DO 2543 I=1,npts SUM=0.d0 DO 2544 J=1,MM IS=J-M-1 II=I+IS IF(II.LT.1)II=iabs(II)+1 IF(II.GT.npts)II=npts-(II-npts-1) SUM=SUM+ftemp(II)*SPOL(J) 2544 CONTINUE freqsc(I)=sum if((i/100000)*100000.eq.i)write(*,'(1H+,i8,$)')i 2543 CONTINUE write(*,'(1H+,i8/)')i c c...Add the correction to the input frequency and output the result c 106 write(*,211) 211 format(1x/' Applying frequency correction to point:'/1x,$) open(2,file=foutsp,status='unknown') do 30 n=1,npts freqcor=freqs(n)+freqsc(n) write(2,32)freqcor,ispec(n),freqs(n) 32 format(f12.5,i8,f14.5) if((n/100000)*100000.eq.n)write(*,'(1H+,i8,$)')n 30 continue close(2) write(*,'(1H+,i8/)')n C C...output calculated spline function for comparison with input points C open(2,file=fsplpts,status='unknown',err=12) c splmax=-1.d20 splmin=1.d20 c nmult=fimult do 20 n=1,npts,nmult write(2,22)freqs(n),freqsc(n) 22 format(f12.3,f12.6) yy=freqsc(n) if(yy.lt.splmin)splmin=yy if(yy.gt.splmax)splmax=yy 20 continue C close(2) C splrang=(splmax-splmin)*1.1d0 splmean=(splmax+splmin)*0.5d0 splmin=splmean-splrang*0.5d0 splmax=splmean+splrang*0.5d0 call wrigle(filgle,fsplpts,fcalpts,finsp,splmin,splmax) <----- c goto 12 c c...save the resulting spectrum in binary c call savesp <----- c 12 stop end C C_____________________________________________________________________________ C C SUBROUTINE SPLINE(fstart) C C DETERMINATION OF CUBIC SPLINE FUNCTION THROUGH A GIVEN SET OF POINTS C C X,Y - vectors containing on input the X and Y coordinates of points C in order of increasing X C A,B,C,D - vectors containing on output coefficients of the cubic spline C components C N - number of points C C Cubic spline component i is only defined in the region X.i to X.i+1, C value of the spline function in this region should be calculated using C the reduced coordinate X-X.i C C FSTART is the normalisation frequency for the X-axis designed C to improve numerical stability C use globdat c implicit integer*4 (I-N) real*8 SU,T,fstart c c...test whether data is uniformly increasing c do 210 i=2,nspl if(x(i).le.x(i-1))then write(*,211)i,x(i)*fstart open(4,file='frecal_dump',status='unknown') do 5000 nnn=1,nspl write(4,5001)x(nnn)*fstart,y(nnn) 5000 continue 5001 format(F15.5,F12.5) close(4) stop endif 210 continue 211 format(1x//1x,60('*')/ * ' ERROR - Frequencies not increasing at point: ',I9/ * ' frequency: ',F15.5 * //9x,'Inspect the FRECAL_DUMP file'/1x,60('*')//) C C Spline function C n=nspl DO 6 I=1,N A(I)=0.d0 B(I)=0.d0 C(I)=0.d0 D(I)=0.d0 6 CONTINUE C DO 1 I=1,N-1 A(I)=X(I+1)-X(I) 1 CONTINUE C DO 2 I=2,N-1 SU=A(I-1)+A(I) B(I)=A(I-1)/SU C(I)=A(I)/SU D(I)=( (Y(I+1)-Y(I))/A(I) - (Y(I)-Y(I-1))/A(I-1) )/SU 2 CONTINUE C B(1)=2.d0 DO 3 I=2,N-1 K=I+1 B(I)=2.d0-B(K)*C(I)/B(I-1) D(K)=D(K)-B(K)*D(I)/B(I-1) 3 CONTINUE C C(N-1)=3.d0*D(N-1)/B(N-2) DO 4 I=1,N-2 K=N-I C(K)=(3.d0*D(K)-C(K)*C(K+1))/B(K-1) 4 CONTINUE C C(1)=0.d0 C(N)=0.d0 DO 5 I=1,N-1 K=I+1 B(I)=(Y(K)-Y(I))/A(I)-A(I)*(C(K)+2.d0*C(I))/3.d0 D(I)=(C(K)-C(I))/3.d0/A(I) A(I)=Y(I) 5 CONTINUE C GRA=B(1) CA=Y(1)-GRA*X(1) T=X(N)-X(N-1) GRB=B(N-1)+2.d0*C(N-1)*T+3.d0*D(N-1)*T*T CB=Y(N)-GRB*X(N) C RETURN END C_____________________________________________________________________________ C C SUBROUTINE VALY(YY,XX,NFUNCT) C C VALY calculates the value of cubic spline function for a given ordinate. C C YY - contains on output the value of the spline function at XX C XX - value of X for evaluation of the spline function C NFUNCT - contains on output the number of spline component used for the C calculation C C in addition: C C A,B,C,D should contain coefficients of the cubic spline components C X should contain the X values at the beginnings of the cubic C components C GRA,CA and GRB,CB should contain the coefficients of the lower and C upper linear end components resp. C use globdat c REAL*8 XX,YY,T C IF(XX.GT.X(1).AND.XX.LT.X(NSPL))GOTO 55 IF(XX.GE.X(NSPL))GOTO 56 YY=XX*GRA+CA NFUNCT=0 GOTO 57 56 YY=XX*GRB+CB NFUNCT=NSPL GOTO 57 C 55 DO 27 I=NSPL-1,1,-1 IF(XX.LT.X(I))GOTO 27 NFUNCT=I GOTO 28 27 CONTINUE C 28 T=XX-X(NFUNCT) YY=A(NFUNCT)+B(NFUNCT)*T+C(NFUNCT)*T*T+D(NFUNCT)*T*T*T C 57 RETURN END C C_____________________________________________________________________________ c subroutine wrigle(filgle,fila,filb,finsp,splmin,splmax) c use specdat c implicit integer*4 (i-n) implicit real*8 (a-h,o-z) c character*100 fila,filb,filgle,finsp character fspan*12,fcent*12 c xcent=(fstart+fend)*0.5d0 xrang=1.1d0*(fend-fstart) c write(fcent,'(f12.1)')xcent 1 if(fcent(1:1).eq.' ')then fcent=fcent(2:len_trim(fcent)) goto 1 endif c write(fspan,'(f12.1)')xrang 2 if(fspan(1:1).eq.' ')then fspan=fspan(2:len_trim(fspan)) goto 2 endif c open(1,file=filgle(1:len_trim(filgle)),status='unknown') c write(1,3)finsp(1:len_trim(finsp)) 3 format('!',80('-')/'! Spline function for spectrum: ',a/ * '!',80('-')/) c write(1,'(a)')'size 29.5 21.0' write(1,'(a)')' ' write(1,'(a)')'lwd=0.04' write(1,'(a)')'set lwidth lwd' write(1,'(a)')'set join round' write(1,'(2a)')'fcent=',fcent write(1,'(2a)')'fspan=',fspan write(1,'(2a)')'fmin=fcent-fspan/2' write(1,'(2a)')'fmax=fcent+fspan/2' c write(1,'(a)')' ' write(1,'(a)')'amove -1.5 1.5' write(1,'(a)')'begin graph' write(1,'(a)')' nobox' write(1,'(a)')' size 35 21' write(1,'(2a)')' xtitle "frequency\,/MHz" hei 0.8 dist 0.4 ', * 'font texcmr' write(1,'(a,a)')' xaxis min fmin max fmax' write(1,'(a)')' xticks length 0.4' write(1,'(a)')' xsubticks length 0.2' write(1,'(a)')' xlabels hei 0.6 ' write(1,'(a)')'! ' write(1,'(a,f14.2,a,f14.2)')' yaxis min ',splmin, * ' max ',splmax write(1,'(a)')' yticks length 0.4' write(1,'(a)')' ysubticks length 0.2' write(1,'(a)')' ylabels hei 0.6' write(1,'(2a)')' ytitle "\Delta\it f\rm\, /MHz" hei 0.8 ', * 'dist .4 font texcmr' write(1,'(a)')'! ' write(1,'(3a)')' data ',fila(1:len_trim(fila)),' d1=c1,c2' write(1,'(a)')' d1 lstyle 1 lwidth lwd color black' c write(1,'(a)')'! ' write(1,'(3a)')' data ',filb(1:len_trim(filb)),' d2=c1,c2' write(1,'(a)')' d2 marker fcircle msize 0.3 color red' write(1,'(a)')' ' c write(1,'(a)')'end graph' c close(1) c return end c C_____________________________________________________________________________ c SUBROUTINE readsp(finp) c c...read the spectrum in standard binary format c use specdat use globdat C integer*2 ita,itb,itc integer*4 n character*100 finp c filnam=finp(1:len_trim(finp)) 1504 OPEN(2,FILE=filnam,FORM='BINARY',STATUS='OLD',ERR=1515) goto 1516 1515 write(*,'(1x/'' Missing or cannot read input spectrum: '',a//)') * finp stop c c...standard binary format c 1516 write(*,'(1x/'' Reading: '',A)')filnam(1:len_trim(filnam)) c READ(2,ERR=503,END=503)COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC, 1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE, 1 SCANSP,PPS,FRMOD,FRAMPL READ(2)ita if(ita.gt.1)then NPTS=ita READ(2)itb,itc ISMALL=itb ILARGE=itc else READ(2)NPTS endif if(iday.lt.1.or.iday.gt.31)goto 503 if(imon.lt.1.or.imon.gt.12)goto 503 if(iyear.lt.1980.or.iyear.gt.2050)goto 503 c do 800 n=1,72 if(ichar(coment(n:n)).lt.32.or.ichar(coment(n:n)).gt.127) * coment(n:n)='#' 800 continue do 801 n=1,20 if(ichar(sample(n:n)).lt.32.or.ichar(sample(n:n)).gt.127) * sample(n:n)='#' 801 continue do 802 n=1,6 if(ichar(lamp(n:n)).lt.32.or.ichar(lamp(n:n)).gt.127) * lamp(n:n)='#' if(ichar(scansp(n:n)).lt.32.or.ichar(scansp(n:n)).gt.127) * scansp(n:n)='#' 802 continue C fstart=0.d0 fend=0.d0 C if(npts.gt.maxpts)then write(*,'(1x//1x,i12,'' points in spectrum, while'', * '' present maximum is'',i8//)')npts,maxpts pause stop endif ISMALL= 2000000000 ILARGE=-2000000000 baslev=0.d0 DO 20 N=1,NPTS READ(2,ERR=494)ita ISPEC(N)=ita if(ispec(n).lt.ismall)ismall=ispec(n) if(ispec(n).gt.ilarge)ilarge=ispec(n) baslev=baslev+ispec(n) 20 CONTINUE imean=(ilarge+ismall)/2 baslev=baslev/real(npts) C goto 495 494 write(*,'(1x//14x, * ''----- PROBLEMS WITH INPUT OF INTENSITIES at pt. '',i5//)')n stop 503 WRITE(*,492) 492 FORMAT(1X/'---- This is not a file in the IFPAN binary format'///) stop c 495 READ(2,err=700,end=700)FSTART,FEND,FINCR WRITE(*,'(1x/'' Frequency from'',F13.2/ * '' to'',F13.2, 1 '', step size'',F10.5)')FSTART,FEND,FINCR IF(FEND.LE.FSTART.or.fincr.lt.0.d0)goto 700 goto 1000 c 700 WRITE(*,'(1x/10x,''PROBLEMS WITH FREQUENCY LIMITS:''//14x, * ''This may be a file in one of the old data standards,'', * /14x,''in which case try using MODSPE''//)') stop C C...fill out the frequency table c 1000 DO 4566 N=1,NPTS FREQS(N)=FSTART+(N-1)*FINCR 4566 CONTINUE c 4565 CLOSE(2) C return end C_____________________________________________________________________________ C C SUBROUTINE SAVESP C C This routine saves the current spectrum into a file: C C a) in FREQLIN type binary format which consists of header, C INTEGER*2 intensities, frequency limits and increment. C b) simple two column ASCII output C C Binary file format byte by byte: C C COMENT character*72 C IDAY,IMON,IYEAR 3 * integer*2 C IHOUR,IMIN,ISEC 3 * integer*2 C LAMP character*6 C VKMST,VKMEND 2 * real*4 C GRID real*4 C SAMPLE character*20 C SAMPRE real*4 C GAIN,TIMEC,PHASE 3 * real*4 C SCANSP character*6 not used in FREQLIN standard C PPS real*4 =PPS+500 in FREQLIN standard C FRMOD,FRAMPL 2 * real*4 C then either: C NPTS integer*2 C ISMALL,ILARGE 2 * integer*2 these are I*4 internally C or C -1 integer*2 C NPTS integer*4 C ISPEC NPTS * integer*2 this is I*4 internally C FSTART,FEND,FINCR 3 * real*8 point f's are determined by FSTART and FINCR C NCALPT integer*2 =0 ensuring backwards compatibility C C USE DFLIB use globdat use specdat C IMPLICIT INTEGER*2 (I-N) IMPLICIT REAL*4 (A-H,O-Z) INTEGER*4 N,ioldr,imean C c REAL*8 rtemp,fout CHARACTER FILOUT*70,CHARFL,STAT*7 C C...preliminaries C 51 write(*,50) 50 format(1x/' The spectrum for rescaling is from file:') write(*,'(5x,a)')filnam(1:len_trim(filnam)) write(*,52) 52 format(1x/' Specify file name for output'/ * ' (no path, ENTER for exit) .... ',$) read(*,'(a)',err=51)filout if(len_trim(filout).lt.1)return C C...Decide on output type and open output file C 31 write(*,30) 30 format(1x/ 5x,'Type of output: 0 = binary'/ * 24x,'1 = two column ASCII'// * 19x,'.... ',$) read(*,'(i5)',err=31)itout if(itout.lt.0.or.itout.gt.1)goto 31 C 22 stat='NEW' 26 if(itout.eq.0)then OPEN(3,FILE=FILOUT(1:len_trim(filout)),STATUS=STAT, * FORM='BINARY',ERR=21) WRITE(3)COMENT,IDAY,IMON,IYEAR,IHOUR,IMIN,ISEC, 1 LAMP,VKMST,VKMEND,GRID,SAMPLE,SAMPRE,GAIN,TIMEC,PHASE, 1 SCANSP,PPS,FRMOD,FRAMPL else OPEN(3,FILE=FILOUT(1:len_trim(filout)),STATUS=STAT,ERR=21) write(3,40)filnam(1:len_trim(filnam)) 40 format('!'/'! Output from spectrum: ',a/'!') endif C c...Intensities of spectral data points - centre intensity on zero c ioldr=ilarge-ismall ISMALL= 1000000000 ILARGE=-1000000000 do 100 n=1,npts if(ispec(n).lt.ismall)ismall=ispec(n) if(ispec(n).gt.ilarge)ilarge=ispec(n) 100 continue imean=(ilarge+ismall)/2 C C...simple ASCII output C if(itout.eq.1)then 34 write(*,33)ismall,ilarge 33 format(1x/ 5x,'Y-axis limits of spectrum =',i12,',',i12// * 5x,' New limits (1/0) ? ',$) read(*,'(i5)',err=4)ichang if(ichang.eq.0)then scale=1.d0 else 37 write(*,36) 36 format(1x/5x,'New limits Ymin,Ymax = ',$) read(*,*,err=37)yminl,ymaxl if(ymaxl.le.yminl)goto 37 scale=(ymaxl-yminl)/real(ilarge-ismall) endif c fout=fstart-fincr do 38 n=1,npts fout=fout+fincr rtemp=yminl+scale*(ispec(n)-ismall) WRITE(3,39)fout,rtemp 38 continue 39 format(f15.5,f20.8) c goto 32 endif C C...standard binary output C 4 write(*,3)ioldr,ilarge-ismall 3 format(1x/ 5x,'Dynamic range of old spectrum =',i12/ * 5x,'Dynamic range of new spectrum =',i12/ * 5x,' Maximum possible = 64000 ', * '(2-byte limit)'// * 5x,' New range (1/0) ? ',$) read(*,'(i5)',err=4)ichang if(ichang.eq.0)then scale=1.d0 else 7 write(*,6) 6 format(1x/5x,' Value of new range = ',$) read(*,*,err=7)rannew scale=rannew/real(ilarge-ismall) endif C if(scale*(ilarge-ismall).gt.64000.) * scale=64000./real(ilarge-ismall) c if(npts.le.32000)then WRITE(3)int2(NPTS),int2(scale*(ISMALL-imean)), * int2(scale*(ILARGE-imean)) else WRITE(3)int2(-1),npts endif c write(*,41) 41 format(1x//' Now working on point:'/1x,$) DO 20 N=1,NPTS WRITE(3)int2(scale*(ispec(n)-imean)) if((n/100000)*100000.eq.n)write(*,'(1H+,i8,$)')n 20 CONTINUE write(*,'(1x)') ncalpt=0 C C...Frequency limits and increment C WRITE(3)FSTART,FEND,FINCR,ncalpt 32 CLOSE(3) write(*,'(1x//5x,''Output has been written to:'')') write(*,'(5x,a)')filout(1:len_trim(filout)) WRITE(*,'(1x/40x,''Press ENTER to continue '',$)') 200 read(*,'(i5)')ik GOTO 23 c 21 WRITE(*,'(1X/'' THIS FILE ALREADY EXISTS, OVERWRITE (Y/N)? '' * ,$)') READ(*,1300,ERR=21)CHARFL IF(CHARFL.EQ.'Y'.OR.CHARFL.EQ.'y')THEN STAT='UNKNOWN' GOTO 26 ENDIF IF(CHARFL.EQ.'N'.OR.CHARFL.EQ.'n')THEN 27 WRITE(*,'('' New file name (or ENTER to exit): '',$)') READ(*,1300,ERR=27)filout if(len_trim(filout).lt.1)return 1300 FORMAT(A) GOTO 22 ENDIF C 23 RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________