	Program LogCL 
c For analysis of logs of cl* experiments (replacement of GNPLT application)
	parameter (ndata=22000)
	implicit real*8 (a-h,o-z)
	dimension BBC(31,ndata),doy(ndata),wynik(10,ndata),Tc(3)
     & ,Trxg(10,ndata)
	character logfile*50,rec*150
     & ,freqc*5,iac*2,pol(ndata)*3
      character cline*127 ! max 127 characters for command line arguments
c	logical istnieje

      call getcl(cline)

c20	idat=22+namel+1		! data start at this record position

	
            if(cline.eq.' ') then
c	logfile='cl10m1tr.log'
c	logfile='\VLBI\10May-J\logs\cl10c2tr.log'
c	logfile='cl10m2-3c128.log'
c	logfile='\VLBI\10May-J\logs\cl10m2tr.log'
c	logfile='cl10m247.log'
c	logfile='\VLBI\10May-J\logs\cl10l2tr.log'
c	logfile='\VLBI\10Sep\cl10l2a.log'
	logfile='cl12k00.log'
		else
	read(cline(1:127),'(a)') logfile
c98	inquire(file=logfile,exist=istnieje,iostat=ios,read=cline)!,err=99
c	print*,cline(1:10),ios,istnieje
		endif

99	open(1,logfile,action='read',status='old')	!,err=98

	nrep=1			! number of repetitions in onoff
	nr=0				! number of data records read
	freqc='bbc01'
	measNo=0

	do 91 iBBClu=1,16	!18! 17 & 18 are if IFa/c
	iBBC=(iBBClu-1)/2+1
	iac(2:2)='l'
	if(mod(iBBClu,2).eq.0) iac(2:2)='u'
		if(iBBC.lt.9) then
	write(freqc(5:5),'(i1)') iBBC
	iac(1:1)=freqc(5:5)	! reading BBC number, 1 to 8
		else
c	write(freqc(5:5),'(a1)') char(2*iBBC+79)	! 'a' or 'c'
	write(freqc(5:5),'(i1)') iBBClu-16
	iac='ia'		! to read IF level (LCP/RCP for BBC01/BBC02)
	if(iBBC.gt.9) iac='ic'
		endif	
	iplace=61+mod(iBBC+1,2)*(82-61)
c2007.153.00:47:10.62:onoff=1,1,75,3,600,ia,ic,1u,2u,3u,4u,5u,6u,7u,8u
c	iplace=61+mod(iBBC-1,4)*(82-61)

30	read(1,'(a)',end=90) rec
c2007.040.16:06:43.55:bbc01=610.99,a,2.0,2.0
		if(rec(22:26).ne.freqc.or.rec(21:21).ne.':') go to 30
c	print*,rec(28:35)
31	read(rec(28:35),'(f6.2,1x,a1)') BBC(16,nr+1),pol(nr+1)(2:2)
ccc	read(rec(37:41),*) BW	
	sigma=dmax1(sigma1,sigma3,sigma4) ! to detect wrong measurements
c	BBC(22,nr)=sigma
c	print*,'OK',nrep,nr
	if(nrep.eq.1) BBC(15,max(nr,1)) = sigma
	sigma1=0
	sigma3=0
	sigma4=0
cc		if(nr.gt.0.and.(BBC(3,nr).le.0.)) then	!or.sigma.gt.100
cc	BBC(16,nr)=BBC(16,nr+1)
cc	pol(nr)(2:2)=pol(nr+1)(2:2)
cc		else
	nr=nr+1
	pol(nr)(3:3)=iac(2:2)
c	if(nr.eq.1) print*,' |'//pol(nr)//'|'
cc		endif
	pol(nr)(1:1)=iac(1:1)
		if(nr.gt.ndata) then
	print*,'  There are more data available. Quitting at:',ndata
	nr=nr-1
	go to 90
		endif
	if(BBC(16,nr).lt.500.) print*,rec
	measNo=0
130	read(1,'(a)',end=90) rec

	if(rec(21:27).eq.':onoff=') read(rec(28:29),*) nrep
c	if(nrep.ne.1) print*,' rep ='
		if(rec(21:30).eq.'#onoff#APR'.and.rec(32:33).eq.iac) then
      read(rec(36:150),'(f8.2,f6.2,f8.1,1x,2f9.5)')
     &  BBC(21,nr),(BBC(id,nr),id=17,20)	 ! center freq is in (21,nr)
	  if(rec(33:33).eq.'l') then
	BBC(16,nr)=-BBC(16,nr) ! negative for LSB (BBC)
	  endif
		endif
c read azimuth & altitude
		if(rec(21:30).eq.'#onoff#VAL') then
	read(rec(43:52),'(2f5.1)') BBC(30,nr),BBC(31,nr) ! Az, El
c	BBC(31,nr) = 90 - BBC(31,nr)	! zenith distance
		endif
	if(idint(BBC(21,nr+1)).eq.8) print*,rec
		if(rec(22:26).eq.freqc.and.rec(21:21).eq.':'
     &.and.(measNo.eq.0.or.measNo.gt.6)) go to 31
		if(rec(21:27).ne.'#onoff#') go to 130
		if(rec(iplace:iplace+1).ne.iac.or.rec(iplace-1:iplace-1).ne.' ')
     & go to 130
	measNo=measNo+1
	if(measNo.gt.15) print*,Nr,measNo
	if(measNo.gt.nrep*5+2) go to 130
	do 80 i=iplace+2,iplace+2+8+9
80	if(rec(i:i).eq.'$') rec(i:i)='6'	! for overflawed readings
c	print*,rec,measNo,nr

	read(rec(iplace+2:iplace+2+8),'(f9.1)') BBC(measNo,nr) !new	 = new
c	print*,iplace,measNo,rec(iplace+2:iplace+2+8), BBC(measNo,nr)
c
c
	read(rec(iplace+11:iplace+19),'(f9.1)') BBC(measNo+22,nr)
	if(measNo.eq.1) read(rec(iplace+11:iplace+19),'(f9.1)') sigma1
	if(measNo.eq.3) read(rec(iplace+11:iplace+19),'(f9.1)') sigma3
	if(measNo.eq.4) read(rec(iplace+11:iplace+19),'(f9.1)') sigma4
	if(iac(1:1).eq.'i') pol(nr)(2:2)=rec(iplace+1:iplace+1)
	if(measNo.gt.7) print*, Measno,' is >7 !!!',nrep,nr
	if(nrep.eq.1)
     &  read(rec(iplace+11:iplace+19),'(f9.1)') BBC(measNo+7,nr) !all sigmas

c2006.135.21:10:13.14
	read(rec,'(i4,1x,i3,2(1x,i2),1x,f5.2)') iyr,idr,ihr,imr,sr
	doy(nr)=idr+((sr/100.+imr)/100.+ihr)/100.
c	if(dabs(BW-.5).gt.0.1d0) doy(nr)=-doy(nr) ! skip other bandwidth data

c	if(nr.eq.1) print*, nr,measNo,freq,pol(nr),BBC(measNo,nr)
c	if(nr.gt.1900) print*,nr, rec(1:50)
c	k=3
c	if(dabs(BBC(k,16)-prev).gt.10d0) print*,nr, prev,BBC(k,16)	! nr.lt.20.and.nr.gt.14)
c	prev=BBC(k,16)
	go to 130	


90	print*,'  No of measurements read in:', nr, measNo,Nbad
     & ,' '//pol(nr),' '//freqc
 	measNo=0
91	rewind(1)
	close(1)

	do 92 i=40,1,-1
	if(logfile(i:i).eq.'.') logfile(i+1:40) ='dat'
	if(logfile(i:i).eq.'\') go to 93 
92	continue
93	if(i.gt.0.and.cline.eq.' ') logfile =logfile(i+1:40) !=' '
	open(2,logfile)


	write(2,'(a,i4,a)') "    DoY",iyr,
     & "       On      On+c     Off+c    "
     & //"Off       Zero     On      On+c    s/Off+c    s/Off     s/On"
     & //"   s/On+c  s/Off+c    s/Off     s/On s3/On+c  BBC  Freq"
     & //"   Trxg  Flux  DPFU    Gain"
	do 100 nrec=1,2! nr
ccc	if(dabs(doy(nrec)).lt..01) go to 100 ! skip wrong meas's in cl07c2
	write(2,'(1x,f13.8,15f9.1,1x,a3,f7.2,f6.2,f6.1,2f8.5,f7.1)')
     & doy(nrec),(BBC(j,nrec),j=1,15),pol(nrec),(BBC(id,nrec),id=16,21)
c	Trxg(1+k,i) =0
100	continue


	do 120 k=0,2,2	! separate off polarizations
	iac(2:2)=char(k+97)	! 'a' or 'c'
	write(2,'(/a,i4,a)')
     & '   CFreq     Tcal    Trxg   TcalK BBC    DoY',iyr,'    Tsys  '
     &//'SEFD   Comp   Tcal(r) Sour   Cal  T3/T1K T2/T2K T7/T3K   '
     &//' S1     S2     C1     C2     C3    dSl   dS+Cl   "0"    '

	Nprn=0
	do 120 i=1,nr
	if(pol(i)(2:2).ne.iac(2:2).and.pol(i)(1:1).ne.iac(2:2))
     & go to 120
	if(dabs(doy(i)).lt.65.05.and.dabs(doy(i)).gt.65.03)go to 120 ! skip wrong meas's in cltest.log
	if(dabs(doy(i)).lt..01) go to 120 ! skip wrong meas's in cl07c2
cOn      On+c     Off+c    Off       Zero     On      On+c     Off+c    Off       On
c 1       2        3        4          5      6         7        8      9         10 
c      On+c    Off+c     Off       On      On+c  BBC  Freq  Tcal   Flux  DPFU    Gain
c       11       12       13       14       15    -     16   17    18     19      20
	On=(BBC(1,i)+BBC(6,i))/(nrep+1)
	Off=BBC(4,i)
	cal = BBC(3,i)
	calo=(BBC(2,i)+BBC(7,i))/(nrep+1)
	calOn_Off= (BBC(2,i)-BBC(1,i) + BBC(3,i)-BBC(4,i)
     & + BBC(7,i)-BBC(6,i))/3
	souOn_Off= (BBC(2,i)-BBC(3,i) + BBC(1,i)-BBC(4,i)
     & + BBC(6,i)-BBC(4,i))/3
		if(nrep.gt.1) then
	On = On+(BBC(10,i)+BBC(14,i))/(nrep+1)
	Off= (Off+BBC(9,i)+BBC(13,i))/nrep
	cal= (cal+BBC(8,i)+BBC(12,i))/nrep
	calo=calo+(BBC(11,i)+BBC(15,i))/(nrep+1)
		endif
	Flux=BBC(18,i)*BBC(20,i)				! Flux density on source
c	fluxU=Flux/max(1d0,(On-Off))				! Jy/count
	fluxU=Flux/max(1d0,(souOn_Off))				! Jy/count
	wynik(4,i)= (calo-On)/max(1d0,(cal - Off))		! Compression (jak w logu!!!)
	wynik(3,i)=(Off-BBC(5,i))*fluxU			! Tsys(Jy) = SEFD
c	wynik(2,i)=(Off-BBC(5,i))*fluxU*BBC(19,i)		! Tsys(Jy)*DPFU
c	wynik(1,i)= (cal - Off)*fluxU*BBC(19,i) 		! Tcal(Jy)*DPFU
	wynik(1,i)= (calOn_Off)*fluxU*BBC(19,i) 		! Tcal(Jy)*DPFU
	wynik(5,i)= (cal - Off)*fluxU*BBC(19,i)/max(1d0,BBC(17,i))! Tcal(K)/Tcal (niemal OK)
	wynik(6,i)=On-Off			! source amplitude
	wynik(7,i)= cal - Off		! cal ampl.
c 1.0303 18.52  254.4 145.330  1.92
	Tc(1)=(BBC(3,i)-BBC(4,i))
     & /max(1d0,(BBC(1,i)-BBC(4,i)))*Flux*BBC(19,i)
		if(nrep.gt.1) then
	Tc(2)=(BBC(8,i)-BBC(9,i))/(BBC(10,i)-BBC(9,i))*Flux*BBC(19,i)
			if(nrep.gt.2)
     & Tc(3)=(BBC(12,i)-BBC(13,i))/(BBC(14,i)-BBC(13,i))*Flux*BBC(19,i)
		else
	 Tc(2)=(BBC(2,i)-BBC(1,i))
     & /max(1d0,(BBC(2,i)-BBC(3,i))*1.05)*Flux*BBC(19,i) ! compression compensation
	 Tc(3)=(BBC(7,i)-BBC(6,i))
     & /max(1d0,(BBC(6,i)-BBC(4,i)))*Flux*BBC(19,i)
		endif

	Tcal=0
	N=0
	do 115 j=1,2 !nrep
	if(Tc(j).gt.30.or.Tc(j).lt.1.) go to 115
	N=N+1
	Tcal=Tcal+Tc(j)
115	continue
	if(N.gt.0) Tcal=Tcal/N
	if(Tcal.gt.0.and.dabs(Tc(3)-Tcal)/Tcal.lt..07d0)
     & Tcal=(Tcal*N+Tc(3))/(N+1)
		if(BBC(22,i).gt.50.or.Tcal.le.0.) then
	Nbad=Nbad+1
	Tcal=0
c	wynik(1,i)=0
		endif
c Not weighted Tcal (same as in GNPLT) ===========
	cal=    BBC(3,i)-BBC(4,i)
	Sour= (BBC(1,i)+BBC(6,i))/2 - BBC(4,i)

	wynik(1,i)=cal/Sour*Flux*BBC(19,i)
	wynik(2,i)=(Off-BBC(5,i))*fluxU*BBC(19,i)		! Tsys(Jy)*DPFU

c  Try to downweight late off- or on-source cases
c
	if((BBC(6,i) - BBC(4,i))/cal.lt..09) then
		 if(BBC(6,i) .lt. BBC(4,i)) then ! third on is off
		BBC(4,i)=BBC(6,i)
		BBC(26,i)=BBC(28,i)
		BBC(28,i)=BBC(28,i)*100
		 else
		BBC(28,i)=BBC(28,i)*10
		 endif
	endif
c
c	BBC(28,i)=BBC(28,i)*10	! wywalam drugie on
	if((BBC(1,i) - BBC(4,i))/cal.lt..02) BBC(23,i)=BBC(23,i)*10
	if(BBC(5,i) .gt. 1d3)  BBC(5,i)=500  !mod(BBC(5,i),1d3) ! wrong zero
	if(Sour/cal.lt..02) BBC(4+22,i)=BBC(4+22,i)*10

c Variance weighted Tcal ===========
	do 110 j=22,22+6
c      if(BBC(j,i).gt.100d0) BBC(j,i)=BBC(j,i)*10  ! Eliminate evidently wrong data
c - not a big deal though
110   BBC(j,i)=BBC(j,i)**2 ! Variances for weights

c see: http://www.astro.uni.torun.pl/~kb/Reports/EVNcal/EVNcalibration.html
	cal=(  (BBC(3,i)-BBC(4,i))/(BBC(3+22,i)+BBC(4+22,i)) +
     & (BBC(2,i)-BBC(1,i))/(BBC(2+22,i)+BBC(1+22,i)) +
     & (BBC(7,i)-BBC(6,i))/(BBC(7+22,i)+BBC(6+22,i))  )/
     & ( 1/(BBC(3+22,i)+BBC(4+22,i)) + 1/(BBC(2+22,i)+BBC(1+22,i))
     & + 1/(BBC(7+22,i)+BBC(6+22,i))  )
	Sour=(  (BBC(1,i)-BBC(4,i))/(BBC(1+22,i)+BBC(4+22,i)) +
     & (BBC(2,i)-BBC(3,i))/(BBC(2+22,i)+BBC(3+22,i)) !*1.05 ! compression compensation
     &+ (BBC(6,i)-BBC(4,i))/(BBC(6+22,i)+BBC(4+22,i))  )/
     & ( 1/(BBC(1+22,i)+BBC(4+22,i)) +1/(BBC(2+22,i)+BBC(3+22,i))
     &+ 1/(BBC(6+22,i)+BBC(4+22,i))  )

c Calculate system temperature from average weighted with variances
      Zero = BBC(5,i)
	var0 = BBC(5+22,i)
	Ts99=(  (BBC(1,i) - Zero - Sour      )/(BBC(1+22,i) + var0)
     &      + (BBC(2,i) - Zero - Sour - cal)/(BBC(2+22,i) + var0)
     &      + (BBC(3,i) - Zero        - cal)/(BBC(3+22,i) + var0)
     &      + (BBC(4,i) - Zero             )*5/(BBC(4+22,i) + var0)
     &      + (BBC(6,i) - Zero - sour      )/(BBC(6+22,i) + var0)
     &      + (BBC(7,i) - Zero - sour - cal)/(BBC(7+22,i) + var0)   )/
     & (  1/(BBC(1+22,i) + var0) + 1/(BBC(2+22,i) + var0) +
     &    1/(BBC(3+22,i) + var0) + 5/(BBC(4+22,i) + var0) +
     &    1/(BBC(6+22,i) + var0) + 1/(BBC(7+22,i) + var0)   )
c Above, 5 at (BBC(4,i) - Zero) is meant to weigh this value as high 
c as the other 5 measurements taken together

c Below, the 1/0.99 coefficient corrects for the fact that Zero
c corresponds to about -20 dB (i.e. 0.01) of Tsys 
	 Tcal=cal/Sour*Flux*BBC(19,i)
	 TsysA=(BBC(4,i) - Zero)/0.99d0/cal*Tcal  ! raw Tsys value 
	 Tsys=Ts99/0.99d0/cal*Tcal 			! weighted Tsys

c=======================================
	 Nprn=Nprn+1
	if(Nprn.lt.3) write(2,'(1x,4f8.2,1x,a3,f12.6, f7.1,f7.2,f7.4,f7.3,
     &2f7.0,3f7.1,20f7.0)') BBC(21,i), Tcal,BBC(17,i),wynik(1,i),pol(i),
     & doy(i),(wynik(j,i),j=2,7), (Tc(j), j=1,3)
     &,BBC(1,i)-BBC(4,i),BBC(6,i)-BBC(4,i),BBC(3,i)-BBC(4,i) ! S1-3
     &,BBC(2,i)-BBC(1,i),BBC(7,i)-BBC(6,i),BBC(6,i)-BBC(1,i) ! C1-3
     &,BBC(7,i)-BBC(2,i),BBC(4,i)-BBC(5,i),BBC(9,i)-BBC(4,i) ! S+C1-3
     &,BBC(22,i)	! max. sigma
	Trxg(1+k,i) = BBC(21,i)
	Trxg(2+k,i) = Tcal
	Trxg(5+k/2,i) = Tsys
	trxg(7+k/2,i) = BBC(31,i)
	Trxg(9+k/2,i) = TsysA
120	continue
150	continue
cOn      On+c     Off+c    Off       Zero     On      On+c     Off+c    Off       On
c 1       2        3        4          5      6         7        8      9         10 
c      On+c    Off+c     Off       On      On+c  BBC  Freq  Tcal   Flux  DPFU    Gain
c       11       12       13       14       15    -     16   17    18     19      20

	iac='lr'
	do 180 k=0,2,2
	write(2,'(/a)') '  f(MHz)  Tcal(K) Pol  Tsys(K) Alt  BBC    TsysA'
	do 180 j=1,nr
	 Fmin=100000
	do 170 i=1,nr
	if(Trxg(1+k,i).gt.0d0.and.Trxg(1+k,i).le.Fmin) then
	 Fmin=Trxg(1+k,i)
	 Nmin=i
	endif
170 	continue
c	if(Trxg(2+k,Nmin).gt.30d0) go to 180 ! pomin Tcal> 30 K
	if(Trxg(2+k,Nmin).gt.100d0) go to 180 ! pomin Tcal> 100 K
	if(Trxg(1+k,Nmin).gt.0d0.and.Trxg(2+k,Nmin).gt.0d0.and. 
     & Trxg(7+k/2,nmin).gt..001d0) write (2,
     &'(f8.2,f9.3,1x,a1,"cp",f9.2,f5.1,a4,f9.2)') Fmin,Trxg(2+k,Nmin),
     & iac(1+k/2:1+k/2),Trxg(5+k/2,nmin),Trxg(7+k/2,nmin),pol(nmin) !	wynik(2,nmin)
     &,Trxg(9+k/2,nmin)
180	Trxg(1+k,Nmin) = 0
	
	print*, Nbad,'  bad measurements detected'
	print*,'         All data stored in the file '//logfile
	close(2)
	end

	real*8 function BBCdB(icount)
	real*8 dB(0:255)
c Table of dB gain setting against gain number for BBCs (256 levels)
c taken from NRAO code subroutine "module/bblvcode.c"
	data dB/
     &-99.99,-36.16,-30.14,-26.62,-24.12,-22.18,-20.60,-19.26,
     &-18.10,-17.08,-16.16,-15.34,-14.58,-13.89,-13.24,-12.64,
     &-12.08,-11.56,-11.06,-10.59,-10.14, -9.72, -9.32, -8.93,
     & -8.56, -8.21, -7.86, -7.54, -7.22, -6.92, -6.62, -6.34,
     & -6.06, -5.79, -5.53, -5.28, -5.04, -4.80, -4.57, -4.34,
     & -4.12, -3.91, -3.70, -3.49, -3.30, -3.10, -2.91, -2.72,
     & -2.54, -2.36, -2.18, -2.01, -1.84, -1.68, -1.52, -1.36,
     & -1.20, -1.05, -0.90, -0.75, -0.60, -0.46, -0.32, -0.18,
     & -0.04,  0.09,  0.23,  0.36,  0.49,  0.61,  0.74,  0.86,
     &  0.98,  1.10,  1.22,  1.34,  1.45,  1.57,  1.68,  1.79,
     &  1.90,  2.01,  2.11,  2.22,  2.32,  2.42,  2.53,  2.63,
     &  2.73,  2.82,  2.92,  3.02,  3.11,  3.21,  3.30,  3.39,
     &  3.48,  3.57,  3.66,  3.75,  3.84,  3.92,  4.01,  4.09,
     &  4.18,  4.26,  4.34,  4.42,  4.50,  4.58,  4.66,  4.74,
     &  4.82,  4.90,  4.97,  5.05,  5.12,  5.20,  5.27,  5.35,
     &  5.42,  5.49,  5.56,  5.63,  5.70,  5.77,  5.84,  5.91,
     &  5.98,  6.05,  6.11,  6.18,  6.25,  6.31,  6.38,  6.44,
     &  6.51,  6.57,  6.63,  6.70,  6.76,  6.82,  6.88,  6.94,
     &  7.00,  7.06,  7.12,  7.18,  7.24,  7.30,  7.36,  7.42,
     &  7.47,  7.53,  7.59,  7.64,  7.70,  7.75,  7.81,  7.86,
     &  7.92,  7.97,  8.03,  8.08,  8.13,  8.19,  8.24,  8.29,
     &  8.34,  8.39,  8.44,  8.50,  8.55,  8.60,  8.65,  8.70,
     &  8.75,  8.80,  8.84,  8.89,  8.94,  8.99,  9.04,  9.08,
     &  9.13,  9.18,  9.23,  9.27,  9.32,  9.37,  9.41,  9.46,
     &  9.50,  9.55,  9.59,  9.64,  9.68,  9.73,  9.77,  9.81,
     &  9.86,  9.90,  9.94,  9.99, 10.03, 10.07, 10.11, 10.16,
     & 10.20, 10.24, 10.28, 10.32, 10.36, 10.40, 10.44, 10.48,
     & 10.52, 10.56, 10.60, 10.64, 10.68, 10.72, 10.76, 10.80,
     & 10.84, 10.88, 10.92, 10.96, 10.99, 11.03, 11.07, 11.11,
     & 11.15, 11.18, 11.22, 11.26, 11.29, 11.33, 11.37, 11.40,
     & 11.44, 11.48, 11.51, 11.55, 11.58, 11.62, 11.65, 11.69,
     & 11.72, 11.76, 11.79, 11.83, 11.86, 11.90, 11.93, 11.97/
	icoun=max0(0,min0(255,icount))
	BBCdB=dB(icoun)
	end
