	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(12,ndata)
	character logfile*50,rec*150 ,freqc*5,iac*2,pol(ndata)*3,pol16*1
      character source*10,cline*127 ! max 127 characters for command line arguments

      call getcl(cline)


	
            if(cline.eq.' ') then
	logfile='cl14l2tr.log'
		else
	read(cline(1:127),'(a)') logfile
		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'		! current BBC name
	measNo=0			! numer pomiaru w onoff



	do 91 iBBClu=1,16		!18! 17 & 18 are if IFa/b
	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 data for BBC number 1 to 8
		else
	write(freqc(5:5),'(i1)') iBBClu-16 ! 1 or 2 for ia i ib
	iac='ia'		! to read IF level (LCP/RCP for BBC01/BBC02)
	if(iBBC.gt.9) iac='ib' ! only for DBBC (otherwise 'ic')
		endif	
	iplace=61+mod(iBBC+1,2)*(82-61) ! kolumna pomiaru np. 1l lub 2l
c2007.153.00:47:10.62:onoff=1,1,75,3,600,ia,ic,1u,2u,3u,4u,5u,6u,7u,8u

30	read(1,'(a)',end=90) rec
c2007.040.16:06:43.55:bbc01=610.99,a,2.0,1
c		if(rec(22:26).ne.freqc.or.rec(21:21).ne.':') go to 30

		if(rec(22:27).eq.freqc//"=") then ! BBC1-4 => rcp
      read(rec(28:35),'(f6.2,1x,a1)') BBC16,pol16 ! BBCfreq,a/c
	pol16='b'
	read(rec(25:26),'(i2)') NrBBC
	if(NrBBC.lt.5) pol16='a'
		endif


	if(rec(21:33).ne.'#onoff#APR '//iac) go to 30
c	read(rec(34:44),*) skyfreq
	backspace(1)
	go to 32

31	sigma=dmax1(sigma1,sigma3,sigma4) ! to detect wrong measurements
	if(nrep.eq.1) BBC(15,max(nr,1)) = sigma
32	sigma1=0
	sigma3=0
	sigma4=0
	nr=nr+1
	pol(nr)(3:3)=iac(2:2)	! l or u sideband
	pol(nr)(2:2)=pol16
	BBC(16,nr)=BBC16	!BBC freq.
	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*," BBCfreq <500: "//rec(1:44),freqc
     &    ,iac, int(BBC(16,nr)),pol(nr)(2:2),nr
	measNo=0

130	read(1,'(a)',end=90) rec
	if(rec(21:27).eq.':'//freqc//"=") then
        read(rec(28:35),'(f6.2,1x,a1)') BBC16,pol16 ! BBCfreq,a/b
	  go to 31
	endif

	if(rec(21:27).eq.':onoff=') read(rec(28:29),*) nrep

	if(nrep.ne.1) print*,' rep ='
		if(rec(21:33).eq.'#onoff#APR '//iac) then
	  if(measNo.eq.7) then
	measNo=0
	nr=nr+1
	pol(nr)(3:3)=iac(2:2)
	pol(nr)(2:2)=pol16
	pol(nr)(1:1)=pol(nr-1)(1:1)
	BBC(16,nr)=BBC16
	  endif
      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 frequency 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
	read(rec(32:41),'(a)') source
c	BBC(31,nr) = 90 - BBC(31,nr)	! zenith distance
		endif

	if(idint(BBC(21,nr+1)).eq.8) print*,rec
	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
c	if(measNo.gt.15) print*,Nr,measNo
c	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

	read(rec(iplace+2:iplace+2+8),'(f9.1)') BBC(measNo,nr)

	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 - time format in logs
	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

	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) ='txt'
	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
	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	!1! separate off polarizations
	iac(2:2)=char(k/2+97)	! 'a' or 'b' !only for DBBC

	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
	read(pol(i)(1:1),'(i1)') NrBBC
c	if(pol(i)(2:2).ne.iac(2:2))	!.and.pol(i)(1:1).ne.iac(2:2) 
	if((NrBBC/5)*2.ne.k)
     & go to 120

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 (nearly 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-source 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	! removing second on-surce
	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 corresponds
c to about -20 dB (i.e. 0.01) of Tsys in analogue BBC (removed for DBBC)
	 Tcal=cal/Sour*Flux*BBC(19,i)
	 TsysA=(BBC(4,i) - Zero)/cal*Tcal   ! /0.99d0 ! raw Tsys value 
	 Tsys=Ts99/cal*Tcal 			! /0.99d0 ! 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
	Trxg(11+k/2,i) = DoY(i)
120	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

c ============= Print out final results =============

	iac='rl'	! tylko dla DBBC (dla analogowych 'lr')

	if(iac.ne.' ') go to 169 ! always skips this printout
c wydruk po jednym BBC obu polaryzacji 
	do 190 k=0,2,2
	write(2,'(/a)') '  f(MHz)  Tcal(K) Pol  Tsys(K) Alt  BBC    TsysA'
	do 190 Nmin=1,nr
	if(pol(nmin)(1:1).ne.'5'.and.pol(nmin)(1:1).ne.'1') go to 190 ! BBC No 1&5
C	if(Trxg(2+k,Nmin).gt.100d0) go to 190 ! pomin Tcal> 100 K
	 Fmin=Trxg(1+k,Nmin)
	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,i5)') 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),Nmin
190	continue


169	write(2,'(/a,i5,a,i5,a)') ' Measurements detected:',nr,',  bad:'
     &,Nbad,'.    Last observed source: '//source

	isort=0 ! 0 - sortowanie wg czestosci, 1 - wg czasu

	do 180 k=0,2,2
	write(2,'(/a,i4)')  
     &'  f(MHz)  Tcal(K) Pol  Tsys(K) Alt  BBC    TsysA   DoY',iyr
	do 180 j=1,nr
	 Fmin=99999.0	! the lowest sky frequency/time first
	do 170 i=1,nr
		if(isort.eq.0) then ! frequency [MHz]
	if(Trxg(1+k,i).gt.0d0.and.Trxg(1+k,i).le.Fmin) then
	 Fmin=Trxg(1+k,i)
	 Nmin=i
	endif
		else		! time [DoY.hhmmss]
	if(Trxg(11+k/2,i).gt.0d0.and.Trxg(11+k/2,i).le.Fmin) then
	 Fmin=Trxg(11+k/2,i)
	 Nmin=i
	endif
		endif
170 	continue
	if(Nminp.eq.Nmin) go to 180
	Nminp=Nmin
	if(Trxg(2+k,Nmin).le.0d0) go to 180 ! pomin Tcal<= 0 K
	if(Trxg(5+k/2,nmin).le.0d0.or.Trxg(5+k/2,nmin).ge.1d4) go to 180 ! pomin 10000K>Tsys<= 0 K
	if(Trxg(2+k,Nmin).ge.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,f11.6)')
     & Trxg(1+k,Nmin),Trxg(2+k,Nmin),iac(1+k/2:1+k/2),	!Fmin
     & Trxg(5+k/2,nmin),Trxg(7+k/2,nmin),pol(nmin),Trxg(9+k/2,nmin)
c     &,Nmin*2+k/2-(k/2)*Nr	 !	wynik(2,nmin)
     &,Trxg(11+k/2,Nmin)	!=Doy(Nmin)	 !	wynik(2,nmin)
 180	Trxg(1+10*isort+k/(1+isort),Nmin) = 0
c 180	Trxg(1+k,Nmin) = 0 or Trxg(11+k/2,Nmin) = 0
	
	print*, Nbad,'  bad measurements detected'
	print*,'         DBBC based results stored in the file '//logfile
	close(2)
	print*,'   Press ENTER to continue'
	read(*,*)
	
	end
