pro norp_rdt2data_2g,rdata,timejo,mvalid,fi,fv $
    ,nocalib=nocalib,ndata=ndata,npa=npa,stat=stt0 $
    ,smax=smax,fimin=fimin,calibration=calibration $
    ,firaw=amp0,fvraw=pol0, org=org


;+
; NAME:
;       NORP_RDT2DATA_2G
;
; PURPOSE:
;       This procedure is for getting 2GHz data from the raw NoRP data
;
; CALLING SEQUENCE:
;       norp_rdt2data_2g,rdata,mvalid,timejo,fi,fv
;
; INPUTS:
;       rdata: structure of the raw data
;       timejo: structure array with 'time internal format'
;               that must be given by 'norp_rdt2timej'.
;
; OPTIONAL INPUT:
;       ndata: 2-elements long-integer array indicating the
;              position of data corresponding to the given timerange
;
; OUTPUTS:
;       mvalid: byte array indicating the data is valid or not.
;              e.g. if the data(123) obtained at time(123)
;                   was the calibration data mvalid(123)=0.
;                   If the data(124) was the real Sun flux mvalid(124)=1.
;       fi:  I (R+L) data array (mfreq,ndata)
;               data(0,*) is the time-sequencial array of the 1st
;               frequency which are read (indicated by 'rdfreq').
;       fv:  V (R-L) data array (mfreq,ndata)
;
; OPTIONAL OUTPUT:
;
; HISTORY:
;       programmed by T. Kosugi in 1984 May as Fortran program
;       for old NoRP data format
;       re-programmed by H.Sekiguchi 1987 July - 1988 Dec.
;
;       1995 July,20    N,Shinox for new format NoRP data
;       95.04.19   ed.  add ATT correction for pol data ,correction data type
;                    changed intarr() to fltarr().,
;       95.09.04 Ver.2.0  M. Nishio and K. Hori
;       95.11.10 Ver.3.0  M. Nishio refernce level at 9.4GHz= amb -> zero
;  <<Caution: ambient level at 9.4GHz is incorrect. Use zero as ref. level>>
;       95.12.12 Ver.3.1  M. Nishio and K. Shibasaki (1GHz ATT 7dB -> 4.3)
;       98.09.03          sekix  & shinox   add 80GHz
;       98.11.20  K. Hori
;       1999-1-12 TY unify program for old/new NoRP data
;       1999-3-12 TY fimin,smax
;       2000-5-09 TY fimin (80->40)
;       2002-5-30 TY revise ATT levels based on Saitoh's memo (2002-05-09)
;       2003-3-12 TY revise ATT levels for new BE
;                 based on Saitoh's memo (2003-02-15)
;	2004-5-11 MS revise ATT levels for new BE
;                 based on Shinohara's comment
;	2021-8-08 MS add mod_fac for correction of time issue in 2021 Jun-July
;-

if not keyword_set(smax) then smax=50
if not keyword_set(fimin) then fimin=40
if not keyword_set(ndata) then begin
  ndata=lonarr(2)
  ndata(0)=0 & ndata(1)=n_elements(timejo)-1
endif

case tag_names(rdata,/structure_name) of
'NORP0' : version =0
'NORP' : version =1
endcase

tamb=35.0
flxref=(tamb+269.0)/5.899      ; SFU of reference device

date=gt_day(timejo(0),/st)

case version of

0: begin   ; old format
print,'Warning: No data in this frequency'
return
end

1: begin   ; new format
    ; ATT levels
    utime0=anytim(date)
    case 1 of
      (utime0 lt anytim('2002-12-04')) : $
             vatt=10.^[0.,0.708,1.404,2.132,2.751,3.453,4.142,4.875]
      (utime0 lt anytim('2003-06-18')) : $
	     vatt=10.^[0.,0.499,0.987,1.49,1.97,2.47,2.95,3.45]
      else : vatt=10.^[0, 0.495,1.011,1.51,2.01,2.51,3.02,3.52]
    endcase

    ;------ The factor for changing the wavelength (MS)
    ;

if not keyword_set(org) then begin
  case 1 of
    (anytim(date) eq anytim('2021-06-12')): mod_fac = 1.004 * 0.94
    (anytim(date) eq anytim('2021-06-13')): mod_fac = 1.008 * 0.94
    (anytim(date) eq anytim('2021-06-14')): mod_fac = 1.012 * 0.94
    (anytim(date) eq anytim('2021-06-15')): mod_fac = 1.016 * 0.94
    (anytim(date) eq anytim('2021-06-16')): mod_fac = 1.02 * 0.94
    (anytim(date) eq anytim('2021-06-17')): mod_fac = 1.024 * 0.94
    (anytim(date) eq anytim('2021-06-18')): mod_fac = 1.028 * 0.94
    (anytim(date) eq anytim('2021-06-19')): mod_fac = 1.032 * 0.94
    (anytim(date) eq anytim('2021-06-20')): mod_fac = 1.036 * 0.94
    (anytim(date) eq anytim('2021-06-21')): mod_fac = 1.04 * 0.94
    (anytim(date) eq anytim('2021-06-22')): mod_fac = 1.044 * 0.94
    (anytim(date) eq anytim('2021-06-23')): mod_fac = 1.048 * 0.94
    (anytim(date) eq anytim('2021-06-24')): mod_fac = 1.052 * 0.94
    (anytim(date) eq anytim('2021-06-25')): mod_fac = 1.056 * 0.94
    (anytim(date) eq anytim('2021-06-26')): mod_fac = 1.06 * 0.94
    (anytim(date) eq anytim('2021-06-27')): mod_fac = 1.064 * 0.94
    (anytim(date) eq anytim('2021-06-28')): mod_fac = 1.068 * 0.94
    (anytim(date) eq anytim('2021-06-29')): mod_fac = 1.072 * 0.94
    (anytim(date) eq anytim('2021-06-30')): mod_fac = 1.077 * 0.94
    (anytim(date) eq anytim('2021-07-01')): mod_fac = 1.083 * 0.94
    (anytim(date) eq anytim('2021-07-02')): mod_fac = 1.085 * 0.94
    (anytim(date) eq anytim('2021-07-03')): mod_fac = 1.088 * 0.94
    (anytim(date) eq anytim('2021-07-04')): mod_fac = 1.093 * 0.94
    (anytim(date) eq anytim('2021-07-05')): mod_fac = 1.097 * 0.94
    (anytim(date) eq anytim('2021-07-06')): mod_fac = 1.101 * 0.94
    (anytim(date) eq anytim('2021-07-07')): mod_fac = 1.105 * 0.94
    (anytim(date) eq anytim('2021-07-08')): mod_fac = 1.109 * 0.94
    
    (anytim(date) ge anytim('2019-05-15')) and $
       (anytim(date) le anytim('2021-06-11')): mod_fac = 0.94  

    (anytim(date) ge anytim('2021-07-09')): mod_fac = 0.94  

;"0.94" is the correction factor for changing the reciver from DSB to SSB 
; 2023/08/07 M.Shimojo 

  else: mod_fac = 1.
  endcase
endif else mod_fac = 1.

iatt=0
npclb=0
dtsky=[30000,0] ; msec
dtamb=[1000,1000] ; msec
dtzero=[1000,1000] ; msec
dtatt=[0,0] ;msec
dtmgn=[600000,1800000] ;msec
flagsky=1
flagamb=2
flagzero=4
end

endcase

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
intego=bytarr(n_elements(rdata))+1b
integ=intego(ndata(0):ndata(1))
timej=timejo(ndata(0):ndata(1))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
amp=rdata.amp2g & pol=rdata.pol2g & stat=rdata.stat2g

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  extract the calibration data

norp_nclb_stat,stat,flagsky,dtsky,timejo,mskyx,nsky,nskyi
norp_nclb_stat,stat,flagamb,dtamb,timejo,mambx,namb,nambi
norp_nclb_stat,stat,flagzero,dtzero,timejo,mzerox,nzero,nzeroi

if keyword_set(nocalib) then begin
  scale=1.
  fisky=0.
  fvsky=0.
  fizero=0.
  fvzero=0.
  goto,L100
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  check the calibration data

if ((mskyx eq 0) or (mambx eq 0) or (mzerox eq 0)) then begin
  print,'***************************************'
  print,'  There is no calibration data 2G ! '
  print,'  # of valid data for SKY ',mskyx
  print,'  # of valid data for AMB ',mambx
  print,'  # of valid data for ZERO ',mzerox
  print,'***************************************'
  ndx=ndata(1)-ndata(0)+1
  fi=fltarr(ndx)
  fv=fltarr(ndx)
  mvalid=bytarr(ndx)
  return
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  integarage the calibration data

norp_clb_integ,mskyx,nsky,npclb,stat,vatt,oksky,intego $
  ,timejo,timejsky_m,amp,fisky_m,pol,fvsky_m,iatt=iatt

norp_clb_integ,mambx,namb,npclb,stat,vatt,okamb,intego $
  ,timejo,timejamb_m,amp,fiamb_m,iatt=iatt

norp_clb_integ,mzerox,nzero,npclb,stat,vatt,okzero,intego $
  ,timejo,timejzero_m,amp,fizero_m,pol,fvzero_m,iatt=iatt

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; check calibration data

whr=where(oksky,moksky)
whr=where(okamb,mokamb)
whr=where(okzero,mokzero)
if ((moksky eq 0) or (mokamb eq 0) or (mokzero eq 0)) then begin
  print,'***************************************'
  print,'  There is no calibration data 2G ! '
  print,'  # of valid data for SKY ',moksky
  print,'  # of valid data for AMB ',mokamb
  print,'  # of valid data for ZERO ',mokzero
  print,'***************************************'
  ndx=ndata(1)-ndata(0)+1
  fi=fltarr(ndx)
  fv=fltarr(ndx)
  mvalid=bytarr(ndx)
  return
endif

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; spline fit

whr=where(oksky,mok)
case 1 of
(mok ge 3) : begin
; fisky=spline(timejsky_m(whr).time,fisky_m(whr),timej.time)
; fvsky=spline(timejsky_m(whr).time,fvsky_m(whr),timej.time)
; res=linfit(timejsky_m(whr).time,fisky_m(whr))
; fisky=res(0)+res(1)*timej.time
; res=linfit(timejsky_m(whr).time,fvsky_m(whr))
; fvsky=res(0)+res(1)*timej.time
  fisky=total(fisky_m(whr))/n_elements(fisky_m(whr))
  fvsky=total(fvsky_m(whr))/n_elements(fvsky_m(whr))
end
(mok eq 2) : begin
  fisky=total(fisky_m(whr))*0.5
  fvsky=total(fvsky_m(whr))*0.5
end
(mok eq 1) : begin
  fisky=(fisky_m(whr))(0)
  fvsky=(fvsky_m(whr))(0)
end
endcase

whr=where(okamb,mok)
case 1 of
(mok ge 3) : begin
; fiamb=spline(timejamb_m(whr).time,fiamb_m(whr),timej.time)
; res=linfit(timejamb_m(whr).time,fiamb_m(whr))
; fiamb=res(0)+res(1)*timej.time
  fiamb=total(fiamb_m(whr))/n_elements(fiamb_m(whr))
end
(mok eq 2) : begin
  fiamb=total(fiamb_m(whr))*0.5
end
(mok eq 1) : begin
  fiamb=(fiamb_m(whr))(0)
end
endcase

whr=where(okzero,mok)
case 1 of
(mok ge 3) : begin
; fizero=spline(timejzero_m(whr).time,fizero_m(whr),timej.time)
; fvzero=spline(timejzero_m(whr).time,fvzero_m(whr),timej.time)
; res=linfit(timejzero_m(whr).time,fizero_m(whr))
; fizero=res(0)+res(1)*timej.time
; res=linfit(timejzero_m(whr).time,fvzero_m(whr))
; fvzero=res(0)+res(1)*timej.time
  fizero=total(fizero_m(whr))/n_elements(fizero_m(whr))
  fvzero=total(fvzero_m(whr))/n_elements(fvzero_m(whr))
end
(mok eq 2) : begin
  fizero=total(fizero_m(whr))*0.5
  fvzero=total(fvzero_m(whr))*0.5
end
(mok eq 1) : begin
  fizero=(fizero_m(whr))(0)
  fvzero=(fvzero_m(whr))(0)
end
endcase

scale=flxref/(fiamb-fisky) ; ??? maybe incorrect

L100:
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;  restore the attainator, and check the valid data

stt0=stat(ndata(0):ndata(1))
amp0=float(amp(ndata(0):ndata(1)))/float(integ)
pol0=float(pol(ndata(0):ndata(1)))/float(integ)

npa=norp_gt_npa(stt0,iatt=iatt)
natti=norp_gt_natti(npa,timej,dtatt)
fi=vatt(npa)*(amp0-fizero)+fizero
fv=vatt(npa)*(pol0-fvzero)+fvzero

fi=scale*(fi-fisky)*mod_fac
fv=scale*(fv-fvsky)*mod_fac

ndx=ndata(1)-ndata(0)+1
mvalid=(norp_mvalid(ndx,natti) $
    and norp_mvalid(ndx,nskyi-ndata(0)) $
    and norp_mvalid(ndx,nambi-ndata(0)) $
    and norp_mvalid(ndx,nzeroi-ndata(0)) $
    and (fi gt fimin) $
    and (scale gt 0) $
    and (scale lt smax))


if 0 then begin
if (ndata(0) eq 0) then begin
whr=min(where(mvalid,cnt))
if cnt ge 1 then begin
tstt=anytim2ints(timej(whr),offset=dtmgn(0)*1.e-3)
mstt=norp_t2i(timej.time,tstt.time)
mvalid(0:mstt)=0b
endif
endif

if (ndata(1) eq n_elements(timejo)-1) then begin
whr=max(where(mvalid,cnt))
if cnt ge 1 then begin
tend=anytim2ints(timej(whr),offset=-dtmgn(1)*1.e-3)
mend=norp_t2i(timej.time,tend.time)
mvalid(mend:*)=0b
endif
endif
endif


calibration={ $
scale:scale $
,oksky:oksky,timejsky_m:timejsky_m,fisky_m:fisky_m,fvsky_m:fvsky_m $
,fisky:fisky,fvsky:fvsky $
,okzero:okzero,timejzero_m:timejzero_m,fizero_m:fizero_m,fvzero_m:fvzero_m $
,fizero:fizero,fvzero:fvzero $
,okamb:okamb,timejamb_m:timejamb_m,fiamb_m:fiamb_m,fiamb:fiamb $
}

return
end