;
;  This IDL procedures convert the native POLDER format into an HDF file
;  Usage : polder2hdf,FILEIN=<filein>,FILEOUT=<fileout.hdf>,ERROR=message
;   where filein is the path to a POLDER standard file, while fileout.hdf is the path to the output HDF file
; If the parameters are not supplied, a dialog asks for the filenames
;
; message is a diagnostic when a problem occured.  0 if succesfull
;
; This software requires the "POLDER define_xxx.pro routines that are distributed with the Anapol IDL software.  One must set
;  IDL path to access these routines.
;
;  Version 1.0   May 26, 2005, FM Breon   fmbreon@cea.fr
;
PRO polder2hdf,FILEIN=filein,FILEOUT=fileout,Error=Error_Message
FORWARD_FUNCTION path_sep,define_product,define_EOS

IF NOT KEYWORD_SET(filein ) THEN filein =DIALOG_PICKFILE(TITLE='POLDER File Selection',FILTER='P*',/MUST_EXIST)
IF (filein EQ '') THEN RETURN
seppos = STRPOS(filein,path_sep(),/REVERSE_SEARCH)

orbit_label = STRMID(filein,seppos+1,15)
filein      = STRMID(filein,0,seppos+16)+'D'     ; "filein" is the datafile (not the POLDER leader file)

IF NOT KEYWORD_SET(fileout) THEN fileout=DIALOG_PICKFILE(TITLE='Select Filename for HDF-EOS file', FILE=orbit_label+'.hdf',/WRITE)
IF (fileout EQ '') THEN RETURN
;
; Opens the data file
;
OPENR,lun,filein,/get_Lun,/SWAP_IF_LITTLE,ERROR=err
IF err NE 0 THEN BEGIN
  Error_Message='The DATA file corresponding to this product was not found'
  RETURN
ENDIF
;
; One verifies that the input file is indeed a standard POLDER product
;
point_lun,lun,8
ref = BYTARR(11) & readu,lun,ref & ref=STRING(ref)
IF (ref NE 'PAST33131CN' AND ref NE 'P2ST33131CN' AND ref NE 'SPG9N122-31') THEN BEGIN
  Error_Message='The file was not recognized as a standard POLDER product'
  RETURN
ENDIF

; Reads the size and number of data records 
point_lun,lun,52
Nrec=0l & length=0l
readu,lun,Nrec,length
;
; Reads filename (into the product) and deduce filetype
point_lun,lun,36
typ_prod = BYTARR(8)  &  readu,lun,typ_prod
typ_prod = STRING(typ_prod)
parasol  = STRMID(typ_prod,1,1) EQ '3'   ;  Recognizes POLDER or Parasol file

IF Parasol THEN print,' Parasol data file' ELSE print,'POLDER data file'
print,' Product type : ',typ_prod
strProd = define_product(typ_prod,length,ERROR=err,PARASOL=parasol)  ; Information about data type
IF err NE 0 THEN BEGIN
  Error_message = 'Error in "Define_Product" procedure'
  RETURN
ENDIF
;
; Reads the equator node.  One centers the projection on this node.
;
IF STRMID(strProd.type,1,1) EQ 3 THEN LonEq = 0. ELSE BEGIN
  fileL = filein
  STRPUT,fileL,'L',seppos+16
  OPENR,lunL,fileL,/GET_LUN,ERROR=err
  IF err NE 0 THEN BEGIN
    Error_Message='The Leader file was not found'
    RETURN
  ENDIF
  point_lun,lunL,180+360+50
  lon = BYTARR(8)  &  readu,lunL,lon
  READS,STRING(lon),LonEq
  free_lun,lunL
ENDELSE
print,'length,Nrec,lonEq:',length,Nrec,lonEq
;
; Reads data
;
data = BYTARR(length,Nrec)
point_lun,lun,180
readu,lun,data
free_lun,lun    ;  Close file
print,' POLDER data is in memory'
;
; Position of pixels on POLDER grid
;
lin    = REFORM(data[ 6,*]*256s + data[ 7,*])
col    = REFORM(data[ 8,*]*256s + data[ 9,*])
lat = 90.-(lin-0.5)/strprod.resol
Nlinmax = 180*StrProd.resol
NcolM = 2*ROUND(Nlinmax*COS(lat*!DTOR))
dist = ROUND(((col-0.5-NlinMax-lonEq*NcolM/360.+1.5*NcolM) MOD NcolM ) -NcolM/2.  -0.5)  ; Dist is the distance from central longitude.

distmin = MIN(dist,MAX=distmax)
Ncol    = distmax-distmin+1
linmin  = MIN(lin,MAX=linmax)
Nlin = linmax-linmin+1
index = (dist-distmin) + LONG(Ncol)*(linmax-lin)
Corners = [distmin,NlinMax/2-linmin-1,distmax+1,NlinMax/2-linmax]   ;  Corners coordinates in pixel units
NcolM=0. & dist=0. & lat=0.                                         ; Saves memory

; These lines to see the map and the box that is used to save the data in HDF file
;WINDOW, XSIZE=800,YSIZE=400
;MAP_SET,0.,lonEq,/SINUSOIDAL,POS=[0.,0.,1.,1.],/NOBORDER,/CONTINENTS
;plot,[-NlinMax,NlinMax],[-NlinMax/2.,NlinMax/2.],/NODATA,XSTYLE=5,YSTYLE=5,/NOERASE
;oplot,Corners[[0,2,2,0,0]],Corners[[3,3,1,1,3]]

;===================================================================
IF FLOAT(!VERSION.RELEASE) LT 5.5 THEN BEGIN
  Error_Message='Version of IDL greater than 5.4 is needed to use the HDF-EOS capabilities'
  RETURN
ENDIF

EOSfileID = EOS_GD_OPEN(fileout,/CREATE)  ; Creates and opens HDF-EOS file
print,'Open HDF-EOS file: ',fileout
IF EOSfileID EQ -1 THEN BEGIN
  Error_Message='Error opening EOS file '+fileout
  RETURN
ENDIF

status=EOS_EH_IDINFO(EOSfileID, HDFfileID, sdID)  

;  HDF_SD_ATTRSET, sdID, 'ATTR_Glob_HDF_SD',  'Ceci est un attribut global HDF_SD user-defined', /STRING

; Defines Grid and projection parameters
strEOS = define_EOS(LonEq,Ncol,Nlin,Corners,strProd.resol,PARASOL=parasol)

; Creates projection grid
gridID = EOS_GD_CREATE(EOSfileID, strEOS.gridname, strEOS.xdim, strEOS.ydim, strEOS.uplft, strEOS.lowrgt)
; Defines projection parameters
status = EOS_GD_DEFPROJ(gridID, strEOS.projname, strEOS.zonecode, strEOS.spherecode, strEOS.projparam)
; Definition du type de pointage sur les pixels.  I AM NOT POSITIVE ABOUT THIS ONE
status = EOS_GD_DEFPIXREG(gridID, strEOS.pixcent)

; Defines dimensions
status = EOS_GD_DEFDIM(gridID, 'XDim', strEOS.xdim)
status = EOS_GD_DEFDIM(gridID, 'YDim', strEOS.ydim)
IF StrProd.Nd_par GT 0 THEN BEGIN
  status = EOS_GD_DEFDIM(gridID, 'ZDim', strEOS.zdim)
  print, 'Dimensions : ', strEOS.xdim, strEOS.ydim, strEOS.zdim
ENDIF ELSE print, 'Dimensions (Ncol,Nlin): ', strEOS.xdim, strEOS.ydim
;==========================================================
; Definition de la liste de parametres NON DIRECTIONNELS
; Ecriture des ATTRIBUTS locaux predefinis
;==========================================================
HDFcode = [99,21,22,24, 5, 6,99,99,99,99,99,99,23,25]  ;  Ce tableau permet de passer des types IDL (de 0 a 15) aux types HDF
idltype = [1,99,99,2,99,99,1,12,99,99,1];  s peut prendre les valeurs -5,-2,1,2, et 5.  ce tableau permet de convertir ces valeurs en type IDL
FOR ip=0, StrProd.Nu_par-1 DO BEGIN
  s = strprod.U_par[ip].size
  IF s GT 5 THEN CONTINUE
  typecode = HDFcode[idltype[s+5]]         ; Format des donnees
  status  = EOS_GD_DEFFIELD   (gridID,(StrProd.U_par[ip]).name, "XDim,YDim", typecode)
  fieldID = HDF_SD_NAMETOINDEX(sdID  ,(StrProd.U_par[ip]).name)  ; Index associe au champ EOS
  sdsID   = HDF_SD_SELECT     (sdID  ,fieldID                 )  ; Identificateur de SDS associe a l'index du champ EOS

; Definitions des attribut locaux "user-defined"
;   HDF_SD_ATTRSET, sdsID, 'NON-DIR Local Attribute', $
;   'This field may be used for user defined local attribute', /STRING

; Ecriture des attributs HDF locaux "predefined"
  cal= DOUBLE((StrProd.U_par[ip]).slope)
  off=-DOUBLE((StrProd.U_par[ip]).offset)/cal
  CalData={    Cal:  cal, $          ;Calibration Factor
           Cal_Err: 0.0D, $          ;Calibration Error
            Offset: off, $           ;Uncalibrated Offset
        Offset_Err: 0.0D, $          ;Uncalibrated Offset Error
          Num_Type: 0L }             ;Number Type of Uncalibrated Data
  HDF_SD_SETINFO, sdsID, FILL=(StrProd.U_par[ip]).dummy, $
                        ;FORMAT=, $
                        ;RANGE=, $
                         LABEL=(StrProd.U_par[ip]).name, $
                         UNIT =(StrProd.U_par[ip]).unit, $
                        ;COORDSYS=, $
                         CALDATA=CalData
  HDF_SD_ENDACCESS, sdsID  ; Fermeture de l'acces au SDS
ENDFOR

;=======================================================
; Definition de la liste de parametres DIRECTIONNELS
; Ecriture des ATTRIBUTS locaux predefinis
;=======================================================
IF StrProd.Nd_par GT 0 THEN BEGIN
  CalData={     Cal: 0.01d, $        ;Calibration Factor
            Cal_Err: 0.0D, $         ;Calibration Error
             Offset: 0.d, $          ;Uncalibrated Offset
         Offset_Err: 0.0D, $         ;Uncalibrated Offset Error
           Num_Type: 0L }            ;Number Type of Uncalibrated Data
  status  = EOS_GD_DEFFIELD(gridID,  'Sol_Zen_Ang', "XDim,YDim,ZDim", 23)
  fieldID = HDF_SD_NAMETOINDEX(sdID,'Sol_Zen_Ang')   ; Index associe au champ EOS
  sdsID   = HDF_SD_SELECT(sdID,fieldID)              ; Identificateur de SDS associe a l'index du champ EOS
  HDF_SD_SETINFO, sdsID, FILL=0l, LABEL='Sol_Zen_Ang',UNIT='deg.', CALDATA=CalData
  status  = EOS_GD_DEFFIELD(gridID,  'View_Zen_Ang', "XDim,YDim,ZDim", 23)
  fieldID = HDF_SD_NAMETOINDEX(sdID,'View_Zen_Ang')  ; Index associe au champ EOS
  sdsID   = HDF_SD_SELECT(sdID,fieldID)              ; Identificateur de SDS associe a l'index du champ EOS  
  HDF_SD_SETINFO, sdsID, FILL=0l, LABEL='View_Zen_Ang',UNIT='deg.', CALDATA=CalData
  status  = EOS_GD_DEFFIELD(gridID,  'Rel_Azim', "XDim,YDim,ZDim", 23)
  fieldID = HDF_SD_NAMETOINDEX(sdID,'Rel_Azim')      ; Index associe au champ EOS
  sdsID   = HDF_SD_SELECT(sdID,fieldID)              ; Identificateur de SDS associe a l'index du champ EOS
  HDF_SD_SETINFO, sdsID, FILL=0l, LABEL='Rel_Azim',UNIT='deg.', CALDATA=CalData
ENDIF
                   
FOR ip=0, StrProd.Nd_par-1 DO BEGIN
  s = strprod.d_par[ip].size
  IF s GT 5 THEN CONTINUE
  typecode = HDFcode[idltype[s+5]]         ; Format des donnees
  status   = EOS_GD_DEFFIELD(gridID,  (strprod.D_par[ip]).name, "XDim,YDim,ZDim", typecode)
  fieldID  = HDF_SD_NAMETOINDEX(sdID, (strprod.D_par[ip]).name)  ; Index associe au champ EOS
  sdsID    = HDF_SD_SELECT(sdID,fieldID)    ; Identificateur de SDS associe a l'index du champ EOS

; Definitions des attribut locaux "user-defined"
;         HDF_SD_ATTRSET, sdsID, 'NON-DIR Local Attribute', $
;         'This field may be used for user defined local attribute', /STRING

; Ecriture des attributs HDF locaux "predefined"
  cal= DOUBLE((StrProd.D_par[ip]).slope)
  off=-DOUBLE((StrProd.D_par[ip]).offset)/cal
  CalData={    Cal: cal, $         ;Calibration Factor
           Cal_Err: 0.0D, $        ;Calibration Error
            Offset: off, $         ;Uncalibrated Offset
        Offset_Err: 0.0D, $        ;Uncalibrated Offset Error
          Num_Type: 0L }           ;Number Type of Uncalibrated Data
  HDF_SD_SETINFO, sdsID, FILL=(StrProd.D_par[ip]).dummy, $
                              ;FORMAT=, $
                              ;RANGE=, $
                               LABEL=(StrProd.D_par[ip]).name, $
                               UNIT =(StrProd.D_par[ip]).unit, $
                              ;COORDSYS=, $
                               CALDATA=CalData
  HDF_SD_ENDACCESS, sdsID  ; Fermeture de l'acces au SDS
ENDFOR
                   
; Operation de detachement/attachement obligatoires
status = EOS_GD_DETACH(gridID)
gridID = EOS_GD_ATTACH(EOSfileID, strEOS.gridname)

;======================================
; Now write non-directional parameters
;======================================
FOR ip=0, StrProd.Nu_par-1 DO BEGIN
 print,'Write ', StrProd.U_par[ip].name
  s = StrProd.U_par[ip].size
  CASE s OF
  1: BEGIN
       sds = BYTARR(Ncol*LONG(Nlin)) + StrProd.U_par[ip].dummy
       sds[index] = data[StrProd.U_par[ip].fi_by,*]
     END
  2: BEGIN
       sds = UINTARR(Ncol*LONG(Nlin)) + StrProd.U_par[ip].dummy
       sds[index] = data[StrProd.U_par[ip].fi_by,*]*256us + data[StrProd.U_par[ip].fi_by+1,*]
     END
 -2: BEGIN
       sds =  INTARR(Ncol*LONG(Nlin)) + StrProd.U_par[ip].dummy
       sds[index] = data[StrProd.U_par[ip].fi_by,*]*256s  + data[StrProd.U_par[ip].fi_by+1,*]
     END
 -5: BEGIN
       sds = BYTARR(Ncol*LONG(Nlin)) + StrProd.U_par[ip].dummy
       sds[index] = data[StrProd.U_par[ip].fi_by,*] MOD 16b
     END
  5: BEGIN
       sds = BYTARR(Ncol*LONG(Nlin)) + StrProd.U_par[ip].dummy
       sds[index] = data[StrProd.U_par[ip].fi_by,*] / 16b
     END
     ELSE : GOTO, jump_loop_u
  ENDCASE
  nsd = WHERE ( sds EQ (StrProd.U_par[ip]).satur OR $
                sds EQ (StrProd.U_par[ip]).unvalid, nb_nsd)
  IF (nb_nsd GT 0) THEN sds[nsd]=(StrProd.U_par[ip]).dummy
  sds = REFORM(sds,Ncol,Nlin)
  status = EOS_GD_WRITEFIELD(gridID, (StrProd.U_par[ip]).name, sds)
  jump_loop_u:
ENDFOR
;======================================
; Now writes directional parameters (if any)
;======================================
IF StrProd.Nd_par GT 0 THEN BEGIN
  print,' Write Solar, view and zenith angle'
  IF parasol THEN Ndir=16 ELSE Ndir=14
  NonDir = REFORM(data[0:StrProd.Nbyte_ND-1,*],StrProd.Nbyte_ND,Nrec)
  data = REFORM(data[StrProd.Nbyte_ND:*,*],StrProd.Nbyte_D,Ndir,Nrec)       ; On extrait les parametres directionnels
  sza = UINTARR(Ndir,Ncol*LONG(Nlin))
  vza = UINTARR(Ndir,Ncol*LONG(Nlin))
  azi = UINTARR(Ndir,Ncol*LONG(Nlin))
  
  CASE StrProd.type OF
    'L1TBG1': BEGIN
                sza[*,index] = (data[ 5,*,*]*256us + data[ 6,*,*])*0.15
                vza[*,index] = (data[ 7,*,*]*256us + data[ 8,*,*])*0.15
                azi[*,index] = (data[ 9,*,*]*256us + data[10,*,*])*0.6
              END
    'L2TLGA': BEGIN
                sza[0,index] = (NonDir[21,*]*256us + NonDir[22,*])*10 & FOR k=1,Ndir-1 DO sza[k,*]=sza[0,*]
                vza[*,index] = (data[ 1,*,*]*256us + data[ 2,*,*])*10
                azi[*,index] = (data[ 3,*,*]*256us + data[ 4,*,*])*10
              END
    'L2TOGA': BEGIN
                sza[0,index] = (NonDir[17,*]*256us + NonDir[18,*])*10 & FOR k=1,Ndir-1 DO sza[k,*]=sza[0,*]
                vza[*,index] = (data[ 0,*,*]*256us + data[ 1,*,*])*10
                azi[*,index] = (data[ 2,*,*]*256us + data[ 3,*,*])*20
              END
    'L2TRGB': BEGIN
                sza[0,index] = ROUND(ACOS(NonDir[20,*]*0.004+0.2)*!RADEG*100.)
                FOR k=1,Ndir-1 DO sza[k,*]=sza[0,*]
                vza[*,index] = data[ 0,*,*]*50
                azi[*,index] = data[ 1,*,*]*150us+18000us
              END
     ELSE:    BEGIN 
                Error_Message='This Directional product cannot be handled.  Contact FM Breon'
                RETURN
              END
  ENDCASE
  sza = REFORM(TRANSPOSE(sza),Ncol,Nlin,Ndir)
  vza = REFORM(TRANSPOSE(sza),Ncol,Nlin,Ndir)
  azi = REFORM(TRANSPOSE(sza),Ncol,Nlin,Ndir)
  status = EOS_GD_WRITEFIELD(gridID, 'Sol_Zen_Ang', sza)
  status = EOS_GD_WRITEFIELD(gridID, 'View_Zen_Ang', vza)
  status = EOS_GD_WRITEFIELD(gridID, 'Rel_Azim', azi)
ENDIF

FOR ip=0, StrProd.Nd_par-1 DO BEGIN
  print,'Write ',StrProd.D_par[ip].name
  s = StrProd.D_par[ip].size
  CASE s OF
  1: BEGIN
       sds = BYTARR(Ndir,Ncol*LONG(Nlin)) + StrProd.D_par[ip].dummy
       sds[*,index] = data[StrProd.D_par[ip].fi_by,*,*]
     END
  2: BEGIN
       sds = UINTARR(Ndir,Ncol*LONG(Nlin)) + StrProd.D_par[ip].dummy
       sds[*,index] = data[StrProd.D_par[ip].fi_by,*,*]*256us + data[StrProd.D_par[ip].fi_by+1,*,*]
     END
 -2: BEGIN
       sds =  INTARR(Ndir,Ncol*LONG(Nlin)) + StrProd.D_par[ip].dummy
       sds[*,index] = data[StrProd.D_par[ip].fi_by,*,*]*256s  + data[StrProd.D_par[ip].fi_by+1,*,*]
     END
 -5: BEGIN
       sds = BYTARR(Ndir,Ncol*LONG(Nlin)) + StrProd.D_par[ip].dummy
       sds[*,index] = data[StrProd.D_par[ip].fi_by,*,*] MOD 16b
     END
  5: BEGIN
       sds = BYTARR(Ndir,Ncol*LONG(Nlin)) + StrProd.D_par[ip].dummy
       sds[*,index] = data[StrProd.D_par[ip].fi_by,*,*] / 16b
     END
     ELSE : GOTO, jump_loop_d
  ENDCASE
  sds = REFORM(TRANSPOSE(sds),Ncol,Nlin,Ndir)
  nsd = WHERE ( sds EQ (StrProd.D_par[ip]).satur OR $
                sds EQ (StrProd.D_par[ip]).unvalid, nb_nsd)
  IF (nb_nsd GT 0) THEN sds[nsd]=(StrProd.D_par[ip]).dummy
  status = EOS_GD_WRITEFIELD(gridID, (StrProd.D_par[ip]).name, sds)
  jump_loop_d:
ENDFOR

;======================================
; Ecriture des attributs GLOBAUX
;======================================
; Ecriture d'attributs globaux EOS "user-defined"
; status = EOS_GD_WRITEATTR(gridID, "ATTR_Glob_EOS_GD", 'Attribut global EOS_GD user-defined', HDF_TYPE=4)
; status = EOS_GD_WRITEATTR(gridID, "SAMPLING", StrProd.wsamp, HDF_TYPE=21)
status = EOS_GD_DETACH(gridID)   ; Fermeture EOS du fichier
status = EOS_GD_CLOSE(EOSfileID)

print,'Succesfull conversion of POLDER file to HDF-EOS'
Error_Message=0
END

;=======================================================================================
; This function returns the HDF-EOS grid and projection parameters
;---------------------------------------------------------------------------------------
FUNCTION define_EOS,LonEq,Ncol,Nlin,Corners,resol,PARASOL=parasol

Radius = 6371007.1809
fac = Radius*2.*!PI/360.

; Definition des param¶eres de grille (sur la zone extraite)
gridname='SINUSOIDAL'
uplft =DBLARR(2)
lowrgt=DBLARR(2)
xdim=Ncol & ydim=Nlin
IF parasol THEN zdim=16 ELSE zdim=14
uplft = DOUBLE(corners[0:1])/resol*fac
lowrgt = DOUBLE(corners[0:1])/resol*fac

; Definition des parametres de projection (tout le globe)
projname=99l
zonecode=0l
spherecode=-1l
pixcent=0l
Sphere = Radius              ; Radius of reference sphere. If zero, 6370997 meters is used
Centmer = LonEq              ; Longitude of the central meridian
FE = 0                       ; False easting in the same units as the semi-major axis
FN = 0                       ; False northing in the same units as the semi-major axis
NZone= 180.*resol            ; Number of equally spaced latitudinal zones (rows)
RFlag = 2                    ; Right jsutify columns flag is used to indicate what to do in
                             ; zones with an odd nombers of columns
projparam = [   Sphere, $
       0, $
       0, $
       0, $
       Centmer, $
       0, $
       FE, $
       FN, $
       NZone, $
       0, $
       RFlag, $
       0, $
       0 ]

; Definition des attributs globaux

StructHDF={gridname:gridname,xdim:xdim,ydim:ydim,zdim:zdim,uplft:uplft,lowrgt:lowrgt, $
           projname:projname,zonecode:zonecode,spherecode:spherecode, $
           projparam:projparam,pixcent:pixcent}

RETURN, StructHDF

END
;========================================================
;  This function returns the directory separator, as a function of the OS system being used
FUNCTION path_sep
;------------------------------------------------------------------------------------------  
CASE !VERSION.OS_FAMILY OF 
    'MacOS'   : sep=':'
    'unix'    : sep='/'
    'Windows' : sep='\'
ENDCASE
RETURN,sep
END