;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;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
;
;Version 1 -- Date: June 18, 2020 
;Initial Version
;
;Version 2 -- June 1, 2026
;1. Fixed bugs in tie-point to full geolocation conversion.
;2. Accomodated changed variable names based on direct broadcast METimage L1B data.

;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_V2_MAIN
  dirIn="C:/Users/wenhui.wang/my.research/data/EPS-SG/db/"
  filename='W_XX-EUMETSAT-Darmstadt,SAT,SGA1-VII-1B-RAD_C_EUMT_20260524173629_G_V_20260524171915_20260524172914_C_N____.nc'
  pathnameIn=dirIn+filename  
  
  timetag=strmid(filename,70,29)
  dirOut=dirIn+'/'+timetag+'_testreader/'
  FILE_MKDIR, dirOut
 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;MetImage constants
  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']
  versionEUMETSAT=2
  rootDataCalibration='data/calibration_data/'
  rootDataMeasurement='data/measurement_data/'
  Zalt=8
  Zact=8
  nDetectors=24 
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
  fileID=H5F_OPEN(pathnameIn)
  
  METIMAGE_READ_DATA_PUB,fileID, 'data/num_chan', temp
  nBands=(size(temp,/dimension))[0]
  METIMAGE_READ_DATA_PUB,fileID, 'data/num_chan_solar', temp
  nSolarBands=(size(temp,/dimension))[0]
  METIMAGE_READ_DATA_PUB,fileID, 'data/num_chan_thermal', temp
  nThermalBands=(size(temp,/dimension))[0]
  
  METIMAGE_READ_DATA_PUB,fileID, 'data/num_pixels', temp
  DIMact=(size(temp,/dimension))[0] 
  
  ;METIMAGE_READ_DATA_PUB,fileID, 'data/num_lines', temp ; should use this, but not working
  METIMAGE_READ_DATA_PUB,fileID, 'data/measurement_data/vii_443', temp 
  DIMalt=(size(temp,/dimension))[1] 
  

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ; read reduced resolution geolocation (in tie point representation)
  ; and convert to full geolocation
  varName=rootDataMeasurement+varNamesGEO[0] ;'latitude'
  METIMAGE_READ_DATA_PUB, fileID, varName, latsReduced
  varName=rootDataMeasurement+varNamesGEO[1] ;'longitude'
  METIMAGE_READ_DATA_PUB, fileID, varName, lonsReduced
  
  geo_grid = interpolate_tiepoints_3d(lonsReduced, latsReduced, DIMact, DIMalt) ;, tie_factor_x, tie_factor_y)
  lons = geo_grid.lons
  lats = geo_grid.lats
  gridALT=geo_grid.gridY
  gridACT=geo_grid.gridX


  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;  
  varName=rootDataMeasurement+varNamesGEO[2] ;'solar_zenith'
  METIMAGE_READ_DATA_PUB, fileID, varName, szasReduced
  szasReducedR=(szasReduced*!PI/180)
  szasR=bilinear(szasReducedR, gridACT,gridALT)  
  cosszas=cos(szasR)
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;   
  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+'band_averaged_solar_irradiance' ;'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 and plot each bands 
  DEVICE, DECOMPOSED=0
  poss = cgLayout([1,1],IXMARGIN=[0,0],OXMARGIN=[0,0],OYMARGIN=[7,0],XGAP=5,YGAP=2)
  ncolors=255
  f=10
  cgDisplay,DIMact/10,DIMalt/10
  pos=poss[*,0]
  
  for indexBandVII=0, nBands-1 DO BEGIN
    bandNameVII=bandNamesVII[indexBandVII]
    varName=rootDataMeasurement+bandNameVII

    if (indexBandVII lt nSolarBands) then begin
      cgLoadCT, 0,ncolors=ncolors
      METIMAGE_READ_DATA_PUB,fileID, varname, radiances
      ;cgImage, radiances

      ;convert radiance to reflectance
      esun=ESUNs[indexBandVII]
      ;For the Metop-SG METImage Level 1B data, the calculation relies on time-dependent,
      ;band-integrated solar irradiance values that have already been adjusted for the exact
      ;orbital position of the Earth at the millisecond the scan was recorded.
      ;refs=!PI*radiances/(cosszas*esun) *(earth_sun_distance_ratio[0])^2
      refs=!PI*radiances/(cosszas*esun)


      listDataOut.add, refs
      
      varNameOut=bandNameVII+'_ref'

      ;plot reflectance
      minvalue=0.0
      maxvalue=1.6
      indices=where(refs lt minvalue,n)
      if (n gt 0) then refs[indices]=minvalue;
      indices=where(refs gt maxvalue,n)
      if (n gt 0) then refs[indices]=maxvalue;
      
      cgImage, refs, minvalue=0, maxvalue=0.8,position=pos,ncolors=ncolors
      posColorBar=[0.1,0.05,0.9,0.075]
      cgCOLORBAR, POSITION=posColorBar,minrange=minvalue, maxrange=maxvalue,title="Reflectance",$
        ncolors=ncolors,charsize=1,ANNOTATECOLOR='black'

    endif else begin
      cgLoadCT, 33, ncolors=ncolors
      ;read radiance and convert radiance to brightness temperatures
      METIMAGE_READ_BT,fileID, varname, radiances,bts
      listDataOut.add, bts
      
      varNameOut=bandNameVII+'_bt'

      meanBT=mean(bts,/NAN)
      sd=stddev(bts,/NAN)
      minvalue=min(bts,/NAN) ;-10 ;-sd
      maxvalue=max(bts,/NAN) ;+10 ;+sd
      
      indices=where(bts lt minvalue,n)
      if (n gt 0) then bts[indices]=minvalue;
      indices=where(bts gt maxvalue,n)
      if (n gt 0) then bts[indices]=maxvalue;

      ;plot BT 
      cgImage, bts, minvalue=minvlaue, maxvalue=maxvalue,position=pos,ncolors=ncolors
      posColorBar=[0.1,0.05,0.9,0.075]
      cgCOLORBAR, POSITION=posColorBar,minrange=minvalue, maxrange=maxvalue,title="Brightness Temperature (K)",$
        ncolors=ncolors,charsize=1,ANNOTATECOLOR='black'

    endelse  

    varNamesOut=[varNamesOut, varNameOut]
    pathImgOut=dirOut+varNameOut+".png"
    write_png, pathImgOut, TVRD(/TRUE)
    print, pathImgOut, " written!"
  ENDFOR
  H5F_CLOSE, FileID

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

FUNCTION interpolate_tiepoints_3d, tie_lons, tie_lats, nx_full, ny_full
  ; =========================================================
  ; STEP 1 CONVERT LATITUDE/LONGITUDE TO Cartesian 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=tie_lats*deg2Rad
  lonsReducedR=tie_lons*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_tie=N*coslat*coslon
  Y_tie=N*coslat*sinlon
  Z_tie=(1-e2)*N*sinlat   ;equation 7 is incorrect

  ; =========================================================
  ; STEP 2: GENERATE THE TARGET INDEX GRIDS (IX, JY)
  ; =========================================================
  size=size(X_tie,/dimensions)
  nx_tie=size[0]
  ny_tie=size[1]

  x_stream = FINDGEN(nx_full) * (FLOAT(nx_tie - 1) / FLOAT(nx_full - 1))
  y_stream = FINDGEN(ny_full) * (FLOAT(ny_tie - 1) / FLOAT(ny_full - 1))

  nZoneALT=fix(ny_full/8)
  nZoneACT=fix(nx_full/8)

  indexStart=0
  offset=0
  for iALT=0, nZoneALT-1 do begin
    indexEnd=indexStart+8
    indices=findgen(9)/8.0+offset
    if (iALT eq nZoneALT-1) then begin
      indexEnd-=1
      indices=indices[0:7]
    endif
    y_stream[indexStart:indexEnd]=indices
    indexStart+=8

    offset++;
    if ((indexStart mod 24) eq 0) then offset++
  endfor

  indexStart=0
  for iACT=0, nZoneACT-1 do begin
    indexEnd=indexStart+8
    indices=findgen(9)/8.0+iACT
    if (iACT eq nZoneACT-1) then begin
      indexEnd-=1
      indices=indices[0:7]
    endif
    x_stream[indexStart:indexEnd]=indices
    indexStart+=8
  endfor

  ; Replicate vectors into 2D meshgrids (IX and JY arrays)
  IX = REBIN(x_stream, nx_full, ny_full, /SAMPLE)
  JY = TRANSPOSE(REBIN(y_stream, ny_full, nx_full, /SAMPLE))

  ; =========================================================
  ; STEP 3: PERFORM 2D INTERPOLATION
  ; =========================================================
  PRINT, 'Interpolating Cartesian components in 3D space...'
  Xfull = BILINEAR(X_tie, IX, JY)
  Yfull = BILINEAR(Y_tie, IX, JY)
  Zfull = BILINEAR(Z_tie, IX, JY)

  ; =========================================================
  ; STEP 4: CONVERT X, Y, Z BACK TO LATITUDE / LONGITUDE
  ; =========================================================
  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


  RETURN, {lons: lonsFull, lats: latsFull,gridX:IX,gridY:JY}
END