;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;MetImage L1B Reader for the version 1 NOAA STAR simulated MetImage L1B files. IDL functions in this file: 
;1. Convert reduced (tie point representation) geolocation to full geolocation
;2. Read radaince and convert to reflectance or brightness temperatures (BT)
;
;Author: Wenhui Wang
;CISESS/ESSIC, University of Maryland - College Park
;5830 University Research Court
;College Park, MD 20740
;
;Date: June 18, 2020
;
;References
;EPS-SG VII Level 1B Product Format Specification
;;references:
;https://books.google.com/books?id=F7jrCAAAQBAJ&pg=PA284&lpg=PA284&dq=non-iterative+approximation+of+Hofmann-Wellenhof&source=bl&ots=zWs5wSrOrF&sig=ACfU3U2cJtcB5k6PQULpCL7EL3FN29Zrfw&hl=en&sa=X&ved=2ahUKEwjx5IbB6NTnAhUhhXIEHRpRDakQ6AEwAnoECAoQAQ#v=onepage&q=non-iterative%20approximation%20of%20Hofmann-Wellenhof&f=false

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


PRO METIMAGE_L1B_READER_MAIN

  ;dirIn='C:\Users\wenhui.wang\my.research\data\EPS-SG\VII\Simulated\UMDSimulated\'
  dirIn='/data/data444/wwang/share/MetImageL1BReader/'
  pathnameIn=dirIn+'W_xx-noaa-star,SAT,SGA1-VII-1B-RAD_C_EUMT_20200222074900_G_D_20200103053500_20200103054000_T_N____.nc'

  dirOut=pathnameIn+'converted\'
  pathnameOut=pathnameIn+'.fullgeo.nc'


  bandNamesVII=['vii_443','vii_555','vii_668','vii_752','vii_763',$
    'vii_865','vii_914','vii_1240','vii_1375','vii_1630',$
    'vii_2250','vii_3740','vii_3959','vii_4050','vii_6725',$
    'vii_7325','vii_8540','vii_10690','vii_12020','vii_13345']
  varNamesGEO=['latitude','longitude', 'solar_zenith','solar_azimuth','sensor_zenith','sensor_azimuth']
 

  nBands=n_elements(bandNamesVII); 20
  nSolarBands=11  ;# of solar bands
  nThermalBands=9   ;# of thermal bands

  ;Initialize grids for geolocation conversion from tie point representation to full geolocation 
  METIMAGE_CALCULATE_RECONSTRUCT_GEOLOCATION_GRIDS, gridALT, gridACT  


  fileID=H5F_OPEN(pathnameIn)
  
  rootDataCalibration='data/calibration_data/'
  rootDataMeasurement='data/measurement_data/'

  ;read reduced resolution geolocation (in tie point representation)
  varName=rootDataMeasurement+varNamesGEO[0] ;'latitude'
  METIMAGE_READ_DATA_PUB, fileID, varName, latsReduced
  varName=rootDataMeasurement+varNamesGEO[1] ;'longitude'
  METIMAGE_READ_DATA_PUB, fileID, varName, lonsReduced
  METIMAGE_RECONSTRUCT_GEOLOCATION, gridALT, gridACT, latsReduced, lonsReduced, lats, lons
  varName=rootDataMeasurement+varNamesGEO[2] ;'solar_zenith'
  METIMAGE_READ_DATA_PUB, fileID, varName, szasReduced
  szas=bilinear(szasReduced, gridALT,gridACT)
  cosszas=cos(szas*!PI/180.0)
   
  listDataOut=list()
  listDataOut.add, lats
  listDataOut.add, lons
  listDataOut.add, szas
  varNamesOut=varNamesGEO[0:2]


  ;read metadata for solar band radiance --> reflectance conversion
  varname=rootDataCalibration+'integrated_solar_irradiance'
  METIMAGE_READ_DATA_PUB,fileID, varname, ESUNs
  varname='status/satellite/earth_sun_distance_ratio'
  METIMAGE_READ_DATA_PUB,fileID, varname, earth_sun_distance_ratio

  ;read solar bands
  for indexBandVII=0, nBands-1 DO BEGIN

    bandNameVII=bandNamesVII[indexBandVII]
    varName=rootDataMeasurement+bandNameVII

    if (indexBandVII lt nSolarBands) then begin
      METIMAGE_READ_DATA_PUB,fileID, varname, radiances
      ;convert radiance to reflectance
      esun=ESUNs[indexBandVII]
      refs=!PI*radiances/(cosszas*esun)*(earth_sun_distance_ratio[0])^2
      listDataOut.add, refs
      varNamesOut=[varNamesOut, bandNameVII+'_ref']
    endif else begin
      ;read radiance and convert radiance to BT
      METIMAGE_READ_BT,fileID, varname, radiances,bts
      listDataOut.add, bts
      varNamesOut=[varNamesOut, bandNameVII+'_bt']
    endelse  
 
  ENDFOR
  H5F_CLOSE, FileID

  ;write geolocation, radiance, ref/BT to a new nc file
  METIMAGE_NETCDF_WRITE_FILE, pathnameOut, varNamesOut,listDataOut

END

;calculate tie point grid indices for full geolocation
PRO METIMAGE_CALCULATE_RECONSTRUCT_GEOLOCATION_GRIDS, gridALT, gridACT
  ;MetImage constants
  Zalt=8  ;pixel per zone in the scan direction
  Zact=8  ;pixel per zone in the track direction
  DIMact=3144  ;# of pixels in scan direction
  DIMalt=4200  ;# of pixels in track direction
  nDetectors=24 ;#of detectors

  nzone_act=DIMact/Zact
  nzone_alt=DIMalt/Zalt
  nScans=DIMalt/nDetectors
  nTiePointsPerScanAlt=nDetectors/Zalt+1

  ;calculate float index value to bilinear interpolation
  gridALT=fltarr(DIMalt,DIMact)
  gridACT=fltarr(DIMalt,DIMact)

  ;alt direction
  array=fltarr(DIMalt)
  fy=indgen(Zalt)*1.0/Zalt
  zone=0
  index1=0
  for iscan=0, nScans-1 do begin
    for j=0, nTiePointsPerScanAlt-2 do begin
      index2=index1+Zalt-1
      array[index1:index2,0]=fy+zone
      zone++
      index1=index2+1
    endfor
    zone++
  endfor
  for i=0, DIMact-1 do begin
    gridALT[*,i]=array
  endfor

  ;act direction
  fx=indgen(DIMact)*1.0/Zact
  for i=0, DIMalt-1 do begin
    gridACT[i,*]=fx
  endfor


END

;convert reduced lats/lons to full lats/lons
PRO METIMAGE_RECONSTRUCT_GEOLOCATION, gridALT, gridACT, latsReduced, lonsReduced, latsFull, lonsFull
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;convert lat lon to X,Y,Z
  deg2Rad=!PI/180
  a=6378137.0
  b=6356752.3

  a2=a*a
  b2=b*b
  e=sqrt(1-b2/a2)
  e2=e*e
  latsReducedR=latsReduced*deg2Rad
  lonsReducedR=lonsReduced*deg2Rad
  sinlat=sin(latsReducedR)
  sinlat2=sinlat*sinlat
  coslat=cos(latsReducedR)
  sinlon=sin(lonsReducedR)
  coslon=cos(lonsReducedR)

  N=a/SQRT(1-e2*sinlat2)
  N1=a/(1-e2*sinlat2)
  ;N=a2/sqrt(a2*coslat*coslat+b2*sinlat*sinlat)
  X=N*coslat*coslon
  Y=N*coslat*sinlon
  Z=(1-e2)*N*sinlat   ;equation 7 is incorrect
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  xFull=bilinear(X, gridALT,gridACT)
  yFull=bilinear(Y, gridALT,gridACT)
  zFull=bilinear(Z, gridALT,gridACT)
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


  xFull2=xFull*xFull
  yFull2=yFull*yFull
  p=sqrt(xFull2+yFull2)
  tmp1=yFull/(xFull+p)
  lonsFullR=2*atan(tmp1)
  lonsFull=lonsFullR/deg2Rad

  e1=sqrt(a2/b2-1) ;equation 9-10
  e12=e1*e1
  theta=atan((zFull*a)/(p*b))
  sintheta=sin(theta)
  sintheta3=sintheta*sintheta*sintheta
  costheta=cos(theta)
  costheta3=costheta*costheta*costheta
  acostheta=acos(theta)
  latsFullR=atan((zFull+e12*b*sintheta3)/(p-e2*a*costheta3))
  latsFull=latsFullR/deg2Rad

END

;read radiances from l1b nc file and convert to brightness temperature
;return both radiances and bts
PRO METIMAGE_READ_BT,fileID, varname, rads,bts
  bandNamesTIR=['vii_3740','vii_3959','vii_4050','vii_6725',$
    'vii_7325','vii_8540','vii_10690','vii_12020','vii_13345']
  bandname=FILE_BASENAME(varname)
  index=where(bandName eq bandNamesTir, n)
  if (n eq 1) then begin
    METIMAGE_READ_DATA_PUB, fileID, varname, rads

    varname='data/calibration_data/channel_cw_thermal'
    METIMAGE_READ_DATA_PUB, fileID, varname, cws
    varname='data/calibration_data/bt_conversion_a'
    METIMAGE_READ_DATA_PUB, fileID, varname, As
    varname='data/calibration_data/bt_conversion_b'
    METIMAGE_READ_DATA_PUB, fileID, varname, Bs

    lamda=cws[index]
    A=As[index]
    B=Bs[index]
    METIMAGE_RADIANCE_2_BT_METHOD1, rads, lamda, A, B, bts
    ;METIMAGE_RADIANCE_2_BT_METHOD2, rads, lamda, bts1
  ENDIF ELSE BEGIN
    bts=!NULL
    print, varname, " not exist!"
  ENDELSE
END

PRO METIMAGE_READ_DATA_PUB, fileID, varname, dataReturned
  dataID=H5D_OPEN(fileID,varName)
  data=H5D_READ(dataID)  ;no need to read

  attrNamesVII=["add_offset","scale_factor","_FillValue","valid_max","valid_min"] ;;;to be cmpleted
  MYH5D_READ_ATTRIBUTES, dataID,attrNamesVII, attrValues
  add_offset=(attrValues[0])[0]
  scale_factor=(attrValues[1])[0]
  fillValue=(attrValues[2])[0]
  validMax=(attrValues[3])[0]
  validMin=(attrValues[4])[0]

  dim=size(data,/dimensions)
  typecode=SIZE(_FillValue,/TYPE)
  if (dim[0] eq 0) then begin
    dim=[1]
  endif
  dataReturned=make_array(dim,type=typecode)

  if (finite(scale_factor)) then begin
    dataReturned[*]=data*scale_factor
  endif else begin
    dataReturned[*]=data
  endelse
  if (finite(add_offset)) then begin
    dataReturned+=add_offset
  endif

  cond=data eq fillValue or data lt validMin or data gt validMax
  indicesBad=where(cond,nBad)
  if (nBad gt 0) then begin
    dataReturned[indicesBad]=!VALUES.F_NAN
  endif


END

PRO MYH5D_READ_ATTRIBUTES, dataID,attrNames, values
  nAvailable=H5A_GET_NUM_ATTRS(dataID)
  namesAvailable=strarr(nAvailable)
  for i=0, nAvailable-1 do begin
    attrid=H5A_OPEN_IDX(dataID, i)
    namesAvailable[i]=H5A_GET_NAME(attrid)
  endfor

  n=n_elements(attrNames)
  values=list()
  for i=0,n-1 do begin
    attrname=attrnames[i]
    index=WHERE(strmatch(namesAvailable,attrname) EQ 1)
    if (index ge 0) then begin
      attrid=H5A_OPEN_NAME(dataID, attrname)
      value=H5A_READ(attrid)
      values.add,value
    endif else begin
      values.add, [!VALUES.F_NAN]
    endelse
  endfor
END




PRO METIMAGE_NETCDF_WRITE_FILE, pathname, varNames,listData

  ndata=size(varNames,/n_elements)

  fid = NCDF_CREATE(pathname, /CLOBBER)
  NCDF_CONTROL, fid, /FILL ; Fill the file with default values.

  dims=size(listData[0],/dimension)
  dimx=dims[0]
  dimy=dims[1]
  ;define dimensions
  xid = NCDF_DIMDEF(fid, 'x', dimx) ; Define the X dimension.
  yid = NCDF_DIMDEF(fid, 'y', dimy) ; Define the y dimension.
  dims=[xid,yid]

  varids=make_array(ndata,/INTEGER)
  for i=0, ndata-1 do begin
    varName=varNames[i]
    data=listData[i]
    varType=TYPENAME(data)

    case varType of
      'BYTE':varids[i]=NCDF_VARDEF(fid, varName, dims, /BYTE)
      'UBYTE':varids[i]=NCDF_VARDEF(fid, varName, dims, /UBYTE)
      'CHAR':varids[i]=NCDF_VARDEF(fid, varName, dims, /CHAR)
      'SHORT':varids[i]=NCDF_VARDEF(fid, varName, dims, /SHORT)
      'USHORT':varids[i]=NCDF_VARDEF(fid, varName, dims, /USHORT)
      'LONG':varids[i]=NCDF_VARDEF(fid, varName, dims, /LONG)
      'ULONG':varids[i]=NCDF_VARDEF(fid, varName, dims, /ULONG)
      'UINT64':varids[i]=NCDF_VARDEF(fid, varName, dims, /UINT64)
      'DOUBLE':varids[i]=NCDF_VARDEF(fid, varName, dims, /DOUBLE)
      'STRING':varids[i]=NCDF_VARDEF(fid, varName, dims, /BYTE)
      ELSE:varids[i]=NCDF_VARDEF(fid, varName, dims, /FLOAT)
    endcase

  endfor
  NCDF_CONTROL, fid, /ENDEF ; Put file in data mode.

  for i=0, ndata-1 do begin
    varid=varids[i]
    data=listData[i]
    NCDF_VARPUT, fid, varid, data
  endfor

  NCDF_CLOSE,fid
END

;Wenhui Wang
;Note:during calculation, radiance and lamda must be converted to m
;ETIMAGE_RADIANCE_2_BT_METHOD1
;Inputs:
;Ls: spectral radiance in W/(m2 sr um)
;Lamda: center wavelength in um
;A, B: correction coefficients bt_conversion_a and bt_consion_b data field
;Output
;bts: brightness temperatures (K)
PRO METIMAGE_RADIANCE_2_BT_METHOD1, Ls, lamda, A, B, bts
  lamda_m=lamda[0]*1.0e-6  ;convert center wavelength from um --> m
  lamda5=lamda_m^5

  c=299792458d
  h=6.626069e-34
  k=1.38065e-23

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  Ls1=Ls*1.0e6  ;where L must be given in W/(m2 sr m) and not in W/(m2 sr μm).

  c1=2.0d*h*c*c
  C2=h*c/k
  tmp=1.0+C1/(Ls1*lamda5)

  xs=C2/(lamda_m*alog(tmp))

  bts=xs*A[0]+B[0]
END