2017-01-09 15:55:42 +00:00
|
|
|
program fcal
|
|
|
|
|
|
|
|
|
|
! Compute Intercept (A) and Slope (B) for a series of FreqCal measurements.
|
|
|
|
|
parameter(NZ=1000)
|
|
|
|
|
implicit real*8 (a-h,o-z)
|
|
|
|
|
real*8 fd(NZ),deltaf(NZ),r(NZ)
|
|
|
|
|
character infile*50
|
|
|
|
|
character line*80
|
|
|
|
|
character cutc*8
|
|
|
|
|
|
|
|
|
|
nargs=iargc()
|
|
|
|
|
if(nargs.ne.1) then
|
|
|
|
|
print*,'Usage: fcal <infile>'
|
|
|
|
|
print*,'Example: fcal fmtave.out'
|
|
|
|
|
go to 999
|
|
|
|
|
endif
|
|
|
|
|
call getarg(1,infile)
|
|
|
|
|
|
|
|
|
|
open(10,file=infile,status='old',err=997)
|
|
|
|
|
open(12,file='fcal.out',status='unknown')
|
|
|
|
|
open(13,file='fcal.plt',status='unknown')
|
|
|
|
|
|
|
|
|
|
i=0
|
|
|
|
|
do j=1,9999
|
|
|
|
|
read(10,1000,end=10) line
|
|
|
|
|
1000 format(a80)
|
|
|
|
|
i0=index(line,' 0 ')
|
|
|
|
|
i1=index(line,' 1 ')
|
|
|
|
|
if(i0.le.0 .and. i1.le.0) then
|
|
|
|
|
read(line,*,err=5) f,df
|
|
|
|
|
ncal=1
|
|
|
|
|
i=i+1
|
|
|
|
|
fd(i)=f
|
|
|
|
|
deltaf(i)=df
|
|
|
|
|
else if(i1.gt.0) then
|
|
|
|
|
i=i+1
|
|
|
|
|
read(line,*,err=5) f,df,ncal,nn,rr,cutc
|
|
|
|
|
fd(i)=f
|
|
|
|
|
deltaf(i)=df
|
|
|
|
|
r(i)=0.d0
|
|
|
|
|
endif
|
|
|
|
|
5 continue
|
|
|
|
|
enddo
|
|
|
|
|
|
|
|
|
|
10 iz=i
|
|
|
|
|
if(iz.lt.2) go to 998
|
|
|
|
|
call fit(fd,deltaf,r,iz,a,b,sigmaa,sigmab,rms)
|
|
|
|
|
|
|
|
|
|
write(*,1002)
|
2017-01-10 14:51:32 +00:00
|
|
|
1002 format(' Freq DF Meas Freq Resid'/ &
|
2017-01-09 15:55:42 +00:00
|
|
|
' (MHz) (Hz) (MHz) (Hz)'/ &
|
2017-01-10 14:51:32 +00:00
|
|
|
'-----------------------------------------')
|
2017-01-09 15:55:42 +00:00
|
|
|
do i=1,iz
|
|
|
|
|
fm=fd(i) + 1.d-6*deltaf(i)
|
|
|
|
|
calfac=1.d0 + 1.d-6*deltaf(i)/fd(i)
|
|
|
|
|
write(*,1010) fd(i),deltaf(i),fm,r(i)
|
|
|
|
|
write(13,1010) fd(i),deltaf(i),fm,r(i)
|
|
|
|
|
1010 format(f8.3,f9.3,f14.9,f9.3,2x,a6)
|
|
|
|
|
enddo
|
|
|
|
|
calfac=1.d0 + 1.d-6*b
|
|
|
|
|
err=1.d-6*sigmab
|
|
|
|
|
|
|
|
|
|
if(iz.ge.3) then
|
|
|
|
|
write(*,1100) a,b,rms
|
|
|
|
|
1100 format(/'A:',f8.2,' Hz B:',f9.4,' ppm StdDev:',f7.3,' Hz')
|
|
|
|
|
if(iz.gt.2) write(*,1110) sigmaa,sigmab
|
|
|
|
|
1110 format('err:',f6.2,9x,f9.4,23x,f13.9)
|
|
|
|
|
else
|
|
|
|
|
write(*,1120) a,b
|
|
|
|
|
1120 format(/'A:',f8.2,' Hz B:',f9.4)
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
write(12,1130) a,b
|
|
|
|
|
1130 format(f10.4)
|
|
|
|
|
|
|
|
|
|
go to 999
|
|
|
|
|
|
|
|
|
|
997 print*,'Cannot open input file: ',infile
|
|
|
|
|
go to 999
|
|
|
|
|
998 print*,'Input file must contain at least 2 valid measurement pairs'
|
|
|
|
|
|
|
|
|
|
999 end program fcal
|
|
|
|
|
|
|
|
|
|
subroutine fit(x,y,r,iz,a,b,sigmaa,sigmab,rms)
|
|
|
|
|
implicit real*8 (a-h,o-z)
|
|
|
|
|
real*8 x(iz),y(iz),r(iz)
|
|
|
|
|
|
|
|
|
|
sx=0.d0
|
|
|
|
|
sy=0.d0
|
|
|
|
|
sxy=0.d0
|
|
|
|
|
sx2=0.d0
|
|
|
|
|
do i=1,iz
|
|
|
|
|
sx=sx + x(i)
|
|
|
|
|
sy=sy + y(i)
|
|
|
|
|
sxy=sxy + x(i)*y(i)
|
|
|
|
|
sx2=sx2 + x(i)*x(i)
|
|
|
|
|
enddo
|
|
|
|
|
delta=iz*sx2 - sx*sx
|
|
|
|
|
a=(sx2*sy - sx*sxy)/delta
|
|
|
|
|
b=(iz*sxy - sx*sy)/delta
|
|
|
|
|
|
|
|
|
|
sq=0.d0
|
|
|
|
|
do i=1,iz
|
|
|
|
|
r(i)=y(i) - (a + b*x(i))
|
|
|
|
|
sq=sq + r(i)**2
|
|
|
|
|
enddo
|
|
|
|
|
rms=0.
|
|
|
|
|
sigmaa=0.
|
|
|
|
|
sigmab=0.
|
|
|
|
|
if(iz.ge.3) then
|
|
|
|
|
rms=sqrt(sq/(iz-2))
|
|
|
|
|
sigmaa=sqrt(rms*rms*sx2/delta)
|
|
|
|
|
sigmab=sqrt(iz*rms*rms/delta)
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
return
|
|
|
|
|
end subroutine fit
|
