; Copyright 2007 IAS

; This file is part of the Planck Sky Model.
;
; The Planck Sky Model is free software; you can redistribute it and/or modify
; it under the terms of the GNU General Public License as published by
; the Free Software Foundation; version 2 of the License.
;
; The Planck Sky Model is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY, without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
; GNU General Public License for more details.
;
; You should have received a copy of the GNU General Public License
; along with the Planck Sky Model. If not, write to the Free Software
; Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

;+
; initializes the galactic part of the PSM by setting the ancillary
; files and parameters
;
;  @param Nside {in}{required}{type=int} Healpix resolution parameter
;  @keyword smallscale{in}{optional}{type=boolean} ...
;  @keyword polarisation{in}{optional}{type=boolean} ...
;  @keyword planckres{in}{optional}{type=boolean} ...
;
;  @history <p> Modified 13/07/04 to include Synchrotron polar fields. </p>
;  @history <p> Modified 08/06/2006 - MAMD : include files in global structure and add smallscale keyword. </p>
;  @history <p> Modified 18/01/2007 - JPB  : include entry for Dobashi Av map used in free_BM.</p>
;  @history <p> Modified 05/03/2007 - MAMD : change dir to fit WG2
;                                     architecture. add smallscale to
;                                     sky_model structure.</p>
;  @history <p> Modified 30/05/2007 - MAMD add polarisation and planckres keywords.</p>
;  @history <p> Modified 28/02/2008 - MAMD update synchrotron polarisation model
;                       (now based on WMAP 23 GHz U and Q).</p>
;  @history <p> MAMD use nside=2048 templates when needed.</p>
;  @history <p> MM add set_ancilfile function to check whether
;  2048 ancillary file are present.</p>
;  @history <p> 25/08/2009: added lmax keyword, Jacques Delabrouille.</p>
;  @history <p> 02/12/2009: former name psm_sys_var_init changed to gal_sys_var_init, Jacques Delabrouille.</p>
;  @history <p> 02/12/2009: changes due to reorganisation of data directory, Jacques Delabrouille.</p>
;-


FUNCTION set_ancilfile , nside, nside_lowres, template_lowres, template_highres

  IF (nside LE nside_lowres) THEN file = template_lowres ELSE file = template_highres
  ok = TEST_PSM_DATA(file, /get)
  IF NOT FILE_TEST(template_highres) THEN BEGIN
     template_highres = template_lowres
     file = template_lowres
     ok = TEST_PSM_DATA(file, /get)
  ENDIF
  IF NOT file_test(file) THEN PRINT, 'WARNING, file '+ file + ' does not exist.'
  RETURN, file
END


PRO gal_sys_var_init,Nside,smallscale=smallscale,polarisation=polarisation,planckres=planckres,lmax=lmax

;-----------------------------------------------
;----------- set Global variables--- -----------
;-----------------------------------------------

; check what is defined
defsysv, '!PSM_DIR', exists=ex_psm_dir
defsysv, '!PSM_TEMPLATE_DIR', exists=ex_psm_template_dir
defsysv, '!PSM_ANCILLARY_DIR', exists=ex_psm_ancillary_dir
defsysv, '!PSMROOT', exists=ex_psmroot

; define PSM code directory
if not keyword_set(ex_psm_dir) and keyword_set(ex_psmroot) then defsysv, '!PSM_DIR', FIX_SEPARATOR(!PSMROOT)+'galactic/'

; define directory for templates (.sav files)
if not keyword_set(ex_psm_template_dir) then begin
   psm_template_dir = getenv('PSM_TEMPLATE_DIR')
   if not keyword_set(psm_template_dir) and keyword_set(ex_psmroot) then defsysv,'!PSM_TEMPLATE_DIR', FIX_SEPARATOR(!PSMROOT)+'data/Intermediate/galactic/' $
   else defsysv,'!PSM_TEMPLATE_DIR', psm_template_dir
endif

; Set filename for ancillary data
if not keyword_set(ex_psm_ancillary_dir) then begin
   dir = getenv('PSM_ANCILLARY_DIR')
   if not keyword_set(dir) and keyword_set(ex_psmroot) then dir = FIX_SEPARATOR(!PSMROOT)+'data/Ancillary' ; WG2 definition
   if not keyword_set(dir) then dir = getenv('TEMPLATE_DIR')                                 ; backward compatibility (obsolete)
   defsysv, '!PSM_ANCILLARY_DIR', dir
endif

;-----------------------------------------------
;----------- set SKY_MODEL structure -----------
;-----------------------------------------------

eptr=ptr_new()  ;Empty pointer

;--------- thermal dust -----------
sfd_st={alpha1:1.5,alpha2:2.6,f1:0.0309,qratio:11.2,T1_av:9.6,T2_av:16.4}  ;FDS model #7
bbm_st={kappas:eptr,freqs:eptr,Ts:eptr}
thermal_dust = {modelnumber:0, $
                sfd_st:sfd_st, $
                SFD_TEMP:[eptr, eptr], $
                SFD_ALPHA:[0.,0.], $
                SFD_NORM:eptr, $
                SFD_NORM_PREFACTOR:[0.,0.], $
                POLAR_FRACTION:[0.,0.], $
                POLAR_GMAP:eptr, $
                POLAR_ANGLE:eptr, $
;--- the following could be removed for memory management
;                I100:eptr, $
;                SFD_RATIO:eptr, $
;                phi_polar:eptr, $
;                delta_gamma_polar:eptr, $
;-----
                bbm_st:bbm_st, $
                TLS_TEMP:eptr}

;-------- spinning dust ------------
dl98={jnu_cnm:eptr,jnu_wnm:eptr, jnu_wim:eptr, $
      jnu_mol:eptr,jnu_drk:eptr, jnu_rn:eptr, jnu_extra:eptr, jnu_freqs:eptr}
spinning_dust = {modelnumber:0, $
                 WCO:eptr, $
                 NHI:eptr, $
                 Halpha:eptr, $
                 Halpha_mask:eptr, $
                 SPIN23GHZ:eptr, $
                 dl98:dl98, $
                 composition:[0.,0.,0.,0.,0.,0.,0.], $
                 SFD_EBV:eptr}
;--------- free-free -------------
free_free = {modelnumber:0, $
             T_e:7000., $         ; electron temperature
             f_d:0.33, $          ; dust fraction (0-1)
             limit:1.0, $         ; limit (mag.) at which to stop correction
             Halpha:eptr, $
             Halpha_mask:eptr, $
             Ebv:eptr, $
             wmap_mem:eptr, $
             mamd2008:eptr}

;--------- synchrotron ------------
synchrotron = {modelnumber:0, $
               T408:eptr, $
               beta:eptr, $
               polQ23GHz:eptr, $
               polU23GHz:eptr, $
               curvature:eptr}
;               gsync:eptr, $
;               polar_sin:eptr, $
;               polar_cos:eptr}

;--------- Ancillary data ----------
;dir = FIX_SEPARATOR(!PSM_ANCILLARY_DIR)
;length = STRLEN(dir)
;dir = STRMID(dir, 0, (length-1)>0)
dir = GETPSMDIR('psmdata')+'ancillary/'

file = SET_ANCILFILE(nside, 1024,dir+'models/gal/SFD/FINK_Rmap_healpix_1024.fits' , dir+'models/gal/SFD/FINK_Rmap_healpix_2048.fits')
SFD_RATIO = {file:file, resolution:7*60., alpha:0.08, beta:1.0, gamma:-3.0, dataptr:eptr}
; cette carte du rapport 240/100 a une resolution variable de ~40' a ~7 degres.

file = set_ancilfile (nside, 1024,dir+'models/gal/SFD/SFD_i100_healpix_1024.fits',dir+'models/gal/SFD/SFD_i100_healpix_2048_2.fits')
SFD_I100 =  {file:file, resolution:9., alpha:0.1, beta:1.0, gamma:-3.0, dataptr:eptr}

file = SET_ANCILFILE(nside, 1024,dir+'models/gal/SFD/SFD_i100_healpix_1024_no_uchii.fits',dir+'models/gal/SFD/SFD_i100_healpix_2048_no_uchii_4.fits')
SFD_I100_NO_HII =  {file:file, resolution:9., alpha:0.1, beta:1.0, gamma:-3.0, dataptr:eptr}

file = SET_ANCILFILE(nside, 256, dir+'models/gal/SFD/onedeg_SFD_Ebv_h256.fits' , dir+'models/gal/SFD/onedeg_SFD_Ebv_h2048.fits')
SFD_EBV =  {file:file, resolution:60., alpha:0.2, beta:1.0, gamma:-3.0, dataptr:eptr}

file = SET_ANCILFILE(nside, 512, dir+'models/gal/MAMD-2009/gmap_thermaldust_3deg-JDprocessed_512.fits', dir+'models/gal/MAMD-2009/gmap_thermaldust_3deg-JDprocessed_2048.fits')
GMAP_POLAR_DUST = {file:file, resolution:3*60., alpha:0.08, beta:0.5, gamma:-2.3, dataptr:eptr}

;file = set_ancilfile (nside, 512, dir+'/SFD/phi_polar_dust_20deg_512.fits', dir+'/SFD/phi_polar_dust_20deg_2048.fits')
;PHI_POLAR_DUST = {file:file, resolution:20*60., alpha:0.0, beta:0.0, gamma:0.0, dataptr:eptr}

file = set_ancilfile (nside, 512, dir+'models/gal/MAMD-2009/delta_gamma_dust_20deg_512_v2.fits', dir+'/Models/ThermalDust/MAMD-2009/delta_gamma_dust_20deg_2048_v2.fits')
DELTA_GAMMA_POLAR_DUST = {file:file, resolution:20*60., alpha:0.0, beta:0.0, gamma:0.0, dataptr:eptr}

file = SET_ANCILFILE(nside, 512, dir+'models/gal/MAMD-2009/spin23_mamd2008_model2-JDprocessed_512.fits', dir+'models/gal/MAMD-2009/spin23_mamd2008_model2-JDprocessed_2048.fits')
SPINNING_23_MAMD2008_MODEL2 = {file:file, resolution:1.1*60., alpha:0.3, beta:1.1, gamma:-3.0, dataptr:eptr}

;file = SET_ANCILFILE(nside, 512, dir+'models/gal/MAMD-2009/spin23_mamd2008_model4-JDprocessed_512.fits', dir+'models/gal/MAMD-2009/spin23_mamd2008_model4-JDprocessed_2048.fits')
;file = SET_ANCILFILE(nside, 512, dir+'models/gal/MAMD-2009/spin23_mamd2008_model4_n512.fits', dir+'models/gal/MAMD-2009/spin23_mamd2008_model4_n2048.fits')
file = SET_ANCILFILE(nside, 512, dir+'models/gal/MAMD-2009/spin23_mamd2008_model4_v2-JDprocessed_512.fits', dir+'models/gal/MAMD-2009/spin23_mamd2008_model4_v2-JDprocessed_2048.fits')
SPINNING_23_MAMD2008_MODEL4 = {file:file, resolution:1.1*60., alpha:0.3, beta:1.1, gamma:-3.0, dataptr:eptr}

file = set_ancilfile (nside, 256, dir+'models/gal/MAMD-2009/sync_beta_model_4_n256_v2_nested.fits', dir+'models/gal/MAMD-2009/sync_beta_model_4_n2048_v2_nested.fits')
MAMD2008_BETA_SYNC_MODEL4 = {file:file, resolution:5*60., alpha:0., beta:0.0, gamma:0.0, dataptr:eptr}

;file = set_ancilfile (nside, 512, dir+'/Synchrotron/beta_sync_mamd_512.fits', dir+'/Synchrotron/beta_sync_mamd_2048.fits')
;MAMD2007_BETA_SYNC = {file:file, resolution:3*60., alpha:0., beta:0.0, gamma:0.0, dataptr:eptr}

file = set_ancilfile (nside, 512, dir+'models/gal/Synchrotron-Curvature/curvature_minus0.3_512.fits', dir+'models/gal/Synchrotron-Curvature/curvature_minus0.3_2048.fits')
SYNCH_CURVATURE_MAP1 = {file:file, resolution:60., alpha:0., beta:0.0, gamma:0.0, dataptr:eptr}

;file = set_ancilfile (nside, 512, dir+'/Synchrotron/gmap_sync23_3deg_512.fits',  dir+'/Synchrotron/gmap_sync23_3deg_2048.fits')
;GSYNC_23GHZ = {file:file, resolution:3*60., alpha:0.17, beta:0.5, gamma:-2.0, dataptr:eptr}

;file = set_ancilfile (nside, 512, dir+'observations/Haslam/lambda_haslam408_dsds_zerolevel_512.fits', dir+'observations/Haslam/lambda_haslam408_dsds_zerolevel_2048.fits')
;LAMBDA_408 = {file:file, resolution:60., alpha:0.05, beta:1.0, gamma:-3.0, dataptr:eptr}
file = set_ancilfile (nside, 512, dir+'observations/Haslam/lambda_haslam408_dsds_zerolevel-JDprocessed_512.fits', dir+'observations/Haslam/lambda_haslam408_dsds_zerolevel-JDprocessed_2048.fits')
LAMBDA_408 = {file:file, resolution:60., alpha:0.05, beta:1.0, gamma:-3.0, dataptr:eptr}

;file = set_ancilfile (nside, 512,dir+'models/gal/MAMD-2009/wmap_polq_23GHz_3deg_512.fits', dir+'models/gal/MAMD-2009/wmap_polq_23GHz_3deg_2048.fits')
file = set_ancilfile (nside, 512,dir+'models/gal/MAMD-2009/wmap7yr_polq_23GHz_3deg_512_v3.fits', dir+'models/gal/MAMD-2009/wmap7yr_polq_23GHz_3deg_2048_v3.fits')
WMAP_POLQ_23GHZ = {file:file, resolution:3*60., alpha:0.8, beta:1.1, gamma:-2.1, dataptr:eptr}
;WMAP_POLQ_23GHZ = {file:file, resolution:3*60., alpha:0.6, beta:1.1, gamma:-2.5, dataptr:eptr}

;file = set_ancilfile (nside, 512, dir+'models/gal/MAMD-2009/wmap_polu_23GHz_3deg_512.fits', dir+'models/gal/MAMD-2009/wmap_polu_23GHz_3deg_2048.fits')
file = set_ancilfile (nside, 512,dir+'models/gal/MAMD-2009/wmap7yr_polu_23GHz_3deg_512_v3.fits', dir+'models/gal/MAMD-2009/wmap7yr_polu_23GHz_3deg_2048_v3.fits')
WMAP_POLU_23GHZ = {file:file, resolution:3*60., alpha:1.0, beta:1.1, gamma:-2.1, dataptr:eptr}
;WMAP_POLU_23GHZ = {file:file, resolution:3*60., alpha:0.6, beta:1.1, gamma:-2.5, dataptr:eptr}

;file = set_ancilfile (nside, 512, dir+'/Synchrotron/Giardino/NewTemplateSFDNs512.fits', dir+'/Synchrotron/Giardino/NewTemplateSFDNs2048.fits')
;GIARDINO_408 = {file:file, resolution:0., alpha:0, beta:0, gamma:0, dataptr:eptr}

;file = set_ancilfile (nside, 512,dir+'/Synchrotron/Giardino/BetaV1N512_MEM_NoFF.fits',dir+'/Synchrotron/Giardino/BetaV1N2048_MEM_NoFF.fits')
;GIARDINO_INDEX = {file:file, resolution:2*60., alpha:0.02, beta:1.0, gamma:-3.0, dataptr:eptr}

;file = set_ancilfile (nside, 512,dir+'/Synchrotron/Giardino/SpectralIndex_Giardino2002.fit',dir+'/Synchrotron/Giardino/SpectralIndex_Giardino2002_2048.fit')
;GIARDINO2002_INDEX = {file:file, resolution:10*60., alpha:0.0, beta:0.0, gamma:0.0, dataptr:eptr}

file = SET_ANCILFILE(nside, 256, dir+'observations/Halpha/onedeg_diff_halpha_256.fits', dir+'observations/Halpha/onedeg_diff_halpha_2048.fits')
;file = set_ancilfile (nside, 512,dir+'/Free-Free/onedeg_diff_halpha.fits', dir+'/Free-Free/onedeg_diff_halpha_2048.fits')
DICKINSON_HALPHA = {file:file, resolution:60., alpha:0.3, beta:1.1, gamma:-2.6, dataptr:eptr}

file = SET_ANCILFILE(nside, 256, dir+'observations/WMAP/MEM/wmap_k_mem_freefree_3yr_v2_256.fits', dir+'observations/WMAP/MEM/wmap_k_mem_freefree_3yr_v2_2048_2.fits')
;file = set_ancilfile (nside, 256, dir+'/Free-Free/wmap_k_mem_freefree_3yr_v2.fits', dir+'/Free-Free/wmap_k_mem_freefree_3yr_v2_2048_2.fits')
FREE_FREE_WMAP_MEM = {file:file, resolution:60., alpha:0.1, beta:1.1, gamma:-2.6, dataptr:eptr}

;file = set_ancilfile (nside, 256, dir+'models/gal/MAMD-2009/free_free_mamd2008.fits', dir+'/Free-Free/free_free_mamd2008.fits')
;FREE_FREE_MAMD2008 = {file:file, resolution:60., alpha:0.1, beta:1.1, gamma:-2.6, dataptr:eptr}

;file = set_ancilfile (nside, 256, dir+'models/gal/MAMD-2009/free_free_mamd2008-JDprocessed_256.fits', dir+'models/gal/MAMD-2009/free_free_mamd2008-JDprocessed_2048.fits')
;FREE_FREE_MAMD2008 = {file:file, resolution:70., alpha:0.1, beta:1.1, gamma:-2.6, dataptr:eptr}
file = set_ancilfile (nside, 256, dir+'models/gal/MAMD-2009/free_free_mamd2008_v2-JDprocessed_256.fits', dir+'models/gal/MAMD-2009/free_free_mamd2008_v2-JDprocessed_2048.fits')
FREE_FREE_MAMD2008 = {file:file, resolution:70., alpha:0.3, beta:1.1, gamma:-2.6, dataptr:eptr}

;file = set_ancilfile (nside, 1024, dir+'/Lambda/H-alpha/Halpha_data_fwhm06_1024.fits',dir+'/Lambda/H-alpha/Halpha_data_fwhm06_2048.fits')
;FINKBEINER_HALPHA = {file:file, resolution:[5., 60.], alpha:0.3, beta:1.0, gamma:-2.6, $
;                     mask:dir+'/Lambda/H-alpha/Halpha_mask_fwhm06_1024.fits', dataptr:eptr}

WCO_DAME = {file:dir+'observations/CO/CO_Dame_data_512.fits', resolution:30., $
            mask:dir+'observations/CO/CO_Dame_mask_512.fits', dataptr:eptr}
; DAME: resolution is variable : 0.125 to 0.5 degree.

NHI_DICKEY_LOCKMAN = {file:dir+'observations/HI/Dickey_Lockman_data_512.fits', $
                      resolution:[30., 60.], $
                      mask:dir+'observations/HI/Dickey_Lockman_mask_512.fits' ,dataptr:eptr}

DOBASHI_AV = {file:dir+'/Extinction/Dobashi/Av_512_NESTED_G_d_MAP.fits', resolution:0, dataptr:eptr}

ANCILLARY = {SFD_RATIO:SFD_RATIO, $
             SFD_I100:SFD_I100, $
             SFD_I100_NO_HII:SFD_I100_NO_HII, $
             SFD_EBV:SFD_EBV, $
             GMAP_POLAR_DUST:GMAP_POLAR_DUST, $
;             PHI_POLAR_DUST:PHI_POLAR_DUST, $
             DELTA_GAMMA_POLAR_DUST:DELTA_GAMMA_POLAR_DUST, $
;             GIARDINO_408:GIARDINO_408, $
;             GIARDINO_INDEX:GIARDINO_INDEX, $
;             GIARDINO2002_INDEX:GIARDINO2002_INDEX, $
;             MAMD2007_BETA_SYNC:MAMD2007_BETA_SYNC, $
             MAMD2008_BETA_SYNC_MODEL4:MAMD2008_BETA_SYNC_MODEL4, $
             SYNCH_CURVATURE_MAP1:SYNCH_CURVATURE_MAP1, $
             LAMBDA_408:LAMBDA_408, $
             WMAP_POLQ_23GHZ:WMAP_POLQ_23GHZ, $
             WMAP_POLU_23GHZ:WMAP_POLU_23GHZ, $
;             GSYNC_23GHZ:GSYNC_23GHZ, $
             DICKINSON_HALPHA:DICKINSON_HALPHA, $
             FREE_FREE_WMAP_MEM:FREE_FREE_WMAP_MEM, $
             FREE_FREE_MAMD2008:FREE_FREE_MAMD2008, $
;             FINKBEINER_HALPHA:FINKBEINER_HALPHA, $
             SPINNING_23_MAMD2008_MODEL2:SPINNING_23_MAMD2008_MODEL2, $
             SPINNING_23_MAMD2008_MODEL4:SPINNING_23_MAMD2008_MODEL4, $
             WCO_DAME:WCO_DAME, $
             NHI_DICKEY_LOCKMAN:NHI_DICKEY_LOCKMAN, $
             DOBASHI_AV:DOBASHI_AV}

; SKY_MODEL STRUCTURE
if not keyword_set(nside) then begin
   nside=128
   print, 'Nside not defined, use default value : ', nside
endif
if not keyword_set(lmax) then begin
   lmax=3*nside
   print, 'Lmax not defined, use default value (3xnside)'
endif

if not keyword_set(smallscale) then smallscale=0
if not keyword_set(polarisation) then polarisation=0
if not keyword_set(planckres) then planckres=5.

sm={thermal_dust:thermal_dust, $
    spinning_dust:spinning_dust, $
    free_free:free_free, $
    synchrotron:synchrotron, $
    Nside:fix(nside), $
    lmax:fix(lmax), $
    smallscale:smallscale, $
    planckresolution:planckres, $
    polarisation:polarisation, $
    ancillary:ancillary}

defsysv,'!sky_model',ptr_new(sm)

heap_gc   ;remove unused pointers

END