modified to cope with spinning dust

Jean-Philippe Bernard
1 parent 98d63e54
Showing 10 changed files with 222 additions and 95 deletions Show diff stats
src/idl/dustem_fit_sed_readme.pro
src/idl/dustem_fit_sed_spinning_readme.pro
src/idl/dustem_read_all_web3p8.pro
src/idl/dustem_read_grain.pro
src/idl/dustem_read_mix.pro
src/idl/dustem_run.pro
src/idl/dustem_set_params.pro
src/idl/dustem_write_all_web3p8.pro
src/idl/dustem_write_grain_web3p8.pro
src/idl/dustem_write_mix.pro
@@ -218,6 +218,11 @@ st=dustem_set_data(sed=spec)
  
 ;== SET THE FITTED PARAMETERS
 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
+;==== below is to reset the defaults keywords for PAH which in newer version of the fortran include spin and mix, ...
+;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
+;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
+(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
+(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
 ;stop
 ;dustem_sort_params
 ;hprint,*(*!dustem_fit).param_descs
@@ -276,15 +281,19 @@ dustem_sed_plot,(*(*!dustem_fit).current_param_values), $
 ;======================================
 ;====You can exit IDL here and re-enter
 ;======================================
-window,1
+window,2
  
 ;=== RESTORE FIT RESULTS
 ;=== initialise dustem
-dustem_init,mode=mode
+dustem_init,mode=use_mode
 ;file=getenv('DUSTEM_RES')+'DUSTEM_fit_example.sav'
 file=!dustem_res+'DUSTEM_fit_example.sav'
 dustem_restore_sed_fit,file
  
+;==== below is to reset the defaults keywords for PAH which in newer version of the fortran include spin and mix, ...
+(*!dustem_params).grains[0].TYPE_KEYWORDS='logn'
+(*!dustem_params).grains[1].TYPE_KEYWORDS='logn'
+
 ;=== recover best fit values
 ;ind=where(*!dustem_fit_chi2 EQ min(*!dustem_fit_chi2))
 ;res=(*!dustem_fit_params);
@@ -195,7 +195,12 @@ st=dustem_set_data(sed=spec)
 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
 ;stop
 ;help,(*!dustem_params).grains[0],/str
-;(*!dustem_params).grains[0].TYPE_KEYWORDS='mix-logn-spin'
+(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-spin'
+(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-spin'
+;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg'
+;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg'
+;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-mix-spin'
+;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-mix-spin'
  
 ;dustem_sort_params
 ;hprint,*(*!dustem_fit).param_descs
@@ -248,6 +253,13 @@ dustem_init,mode=mode
 file=!dustem_res+'DUSTEM_fit_example.sav'
 dustem_restore_sed_fit,file
  
+;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-chrg'
+;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-chrg'
+;(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-mix-spin'
+;(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-mix-spin'
+(*!dustem_params).grains[0].TYPE_KEYWORDS='logn-spin'
+(*!dustem_params).grains[1].TYPE_KEYWORDS='logn-spin'
+
 ;=== recover best fit values
 ;ind=where(*!dustem_fit_chi2 EQ min(*!dustem_fit_chi2))
 ;res=(*!dustem_fit_params);
-FUNCTION dustem_read_all_web3p8,dir_in,silent=silent
+FUNCTION dustem_read_all_web3p8,dir_in,silent=silent,help=help
+
+;+
+; NAME:
+;    dustem_read_all_web3p8
+; PURPOSE:
+;    Reads all Dustem input data, as of the web available dustem fortran version
+; CATEGORY:
+;    Dustem
+; CALLING SEQUENCE:
+;    st=dustem_read_all_web3p8(dir_in[,/silent])
+; INPUTS:
+;    dir_in: input directory
+; OPTIONAL INPUT PARAMETERS:
+;    None
+; OUTPUTS:
+;    None
+; OPTIONAL OUTPUT PARAMETERS:
+;    None
+; ACCEPTED KEY-WORDS:
+;    help      = If set, print this help
+;    sielnt    = If set, keep silent
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    The dustem idl wrapper must be installed
+; PROCEDURE:
+;    None
+; EXAMPLES
+;    
+; MODIFICATION HISTORY:
+;    Written by J.-Ph. Bernard
+;    see evolution details on the dustem cvs maintained at CESR
+;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
+;-
+
+IF keyword_set(help) THEN BEGIN
+  doc_library,'dustem_read_all_web3p8'
+  st=0
+  goto,the_end
+ENDIF
  
 dir_in_dat=dir_in+'/data/'
 dir_in_qabs=dir_in+'/oprop/'
@@ -13,6 +55,9 @@ Ngrains=n_elements(st_grains)
 file=dir_in_dat+'ISRF.DAT'
 st_isrf=dustem_read_isrf(file,silent=silent,Nisrf=Nisrf)
  
+file=dir_in_dat+'GAS.DAT'
+st_gas=dustem_read_gas(file,silent=silent)
+
 ;=== Read Qabs
  
 st_qabs=ptrarr(Ngrains)
@@ -46,6 +91,7 @@ FOR i=0L,Ngrains-1 DO BEGIN
   ENDELSE
 ENDFOR
  
+;stop
 ;=== Read MIX files
 st_mix = ptrarr(Ngrains)
 FOR i=0L,Ngrains-1 DO BEGIN
@@ -58,6 +104,30 @@ FOR i=0L,Ngrains-1 DO BEGIN
    ENDELSE
 ENDFOR
  
+;=== Read CHRG files
+st_chrg = ptrarr(Ngrains)
+FOR i=0L,Ngrains-1 DO BEGIN
+   IF stregex(st_grains(i).type_keywords, 'chrg', /bool) THEN BEGIN
+      chrg_file=dir_in_dat+'CHRG_'+st_grains(i).grain_type+'.DAT'
+      st=dustem_read_chrg(chrg_file,silent=silent)
+      st_chrg(i)=ptr_new(st)
+   ENDIF ELSE BEGIN
+      st_chrg(i)=ptr_new()
+   ENDELSE
+ENDFOR
+
+;=== Read SPIN files
+st_spin = ptrarr(Ngrains)
+FOR i=0L,Ngrains-1 DO BEGIN
+   IF stregex(st_grains(i).type_keywords, 'spin', /bool) THEN BEGIN
+      chrg_file=dir_in_dat+'SPIN_'+st_grains[i].grain_type+'.DAT'
+      st=dustem_read_spin(chrg_file,silent=silent)
+      st_spin[i]=ptr_new(st)
+   ENDIF ELSE BEGIN
+      st_spin[i]=ptr_new()
+   ENDELSE
+ENDFOR
+
 ;stop
  
 ;=== Read POL and Qabs files
@@ -121,17 +191,11 @@ if !run_tls then begin
    ENDFOR
 ENDIF
  
-;VG Add NH
-;if !run_pol then begin
-;	st={grain:st_grains,isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,align:st_align,qpol:st_qpol,qcirc:st_qcirc,qH:st_qH,tls:st_tls}
-;endif else begin
-;	st={grain:st_grains,isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,tls:st_tls}
-;endelse
-
-;st={Ngrains:Ngrains,G0:G0,Keywords:key_str,grains:st_grains, $
-;    isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,pol:st_pol,qabspol:st_qabspol,tls:st_tls}
+;=== This will become !dustem_params
 st={Ngrains:Ngrains,G0:G0,Keywords:key_str,grains:st_grains, $
-    isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,pol:st_pol,tls:st_tls}
+    isrf:st_isrf,qabs:st_qabs,calor:st_calor,lambda:st_lambda,size:st_size,mix:st_mix,chrg:st_chrg,gas:st_gas,spin:st_spin,pol:st_pol,tls:st_tls}
+
+the_end:
  
 RETURN,st
  
-FUNCTION dustem_read_grain,file,silent=silent,key_str=key_str,G0=G0
+FUNCTION dustem_read_grain,file,silent=silent,key_str=key_str,G0=G0,help=help
+
+;+
+; NAME:
+;    dustem_read_grain
+; PURPOSE:
+;    reads the file GRAIN.dat and returns the corresponding grain structure
+; CATEGORY:
+;    Dustem
+; CALLING SEQUENCE:
+;    st=dustem_read_grain(file,[/silent][,key_str=][,G0=][,/help])
+; INPUTS:
+;    file = file name to be read
+; OPTIONAL INPUT PARAMETERS:
+;    None
+; OUTPUTS:
+;    st: output structure
+; OPTIONAL OUTPUT PARAMETERS:
+;    key_str : a string containing grain keywords
+;    G0      : The G0 value (scaling factor for radiation field)
+; ACCEPTED KEY-WORDS:
+;    help      = If set, print this help
+;    silent      = If set, keeps quiet
+; COMMON BLOCKS:
+;    None
+; SIDE EFFECTS:
+;    None
+; RESTRICTIONS:
+;    The dustem idl wrapper must be installed
+; PROCEDURE:
+;    None
+; EXAMPLES
+;    
+; MODIFICATION HISTORY:
+;    Written by J.-Ph. Bernard, N. Flagey, D. Paradis Jan 2007
+;    see evolution details on the dustem cvs maintained at CESR
+;    Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
+;-
  
 ;=== Format for DUSTEM_WHICH=VERSTRAETE
 ;; # DUSTEM: definition of grain populations
@@ -34,16 +71,19 @@ FUNCTION dustem_read_grain,file,silent=silent,key_str=key_str,G0=G0
 ;; 1.000000
 ;; PAH0              10 mix-logn           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
 ;; PAH1              10 mix-logn           7.8000E-04  2.2400E+00  3.5000E-08  1.2000E-07  6.4000E-08  1.0000E-01
-;; amCBEx            15 logn         	   1.6500E-04  1.8100E+00  6.0000E-08  2.0000E-06  2.0000E-07  3.5000E-01
-;; amCBEx  	     25 plaw-ed		   1.4500E-03  1.8100E+00  4.0000E-07  2.0000E-04 -2.8000E+00  1.5000E-05  1.5000E-05  2.0000E+00
-;; aSil              25 plaw-ed   	   7.8000E-03  3.5000E+00  4.0000E-07  2.0000E-04 -3.4000E+00  2.0000E-05  2.0000E-05  2.0000E+00
+;; amCBEx            15 logn             1.6500E-04  1.8100E+00  6.0000E-08  2.0000E-06  2.0000E-07  3.5000E-01
+;; amCBEx            25 plaw-ed      1.4500E-03  1.8100E+00  4.0000E-07  2.0000E-04 -2.8000E+00  1.5000E-05  1.5000E-05  2.0000E+00
+;; aSil              25 plaw-ed        7.8000E-03  3.5000E+00  4.0000E-07  2.0000E-04 -3.4000E+00  2.0000E-05  2.0000E-05  2.0000E+00
+
+IF keyword_set(help) THEN BEGIN
+  doc_library,'dustem_read_grain'
+  full_st=0.
+  goto,the_end
+ENDIF
  
-;VG
 defsysv, '!run_pol', 0.         ; Global flag to know if polarisation is being run
-;DP
-defsysv,'!run_tls',0.
+defsysv,'!run_tls',0.           ; Global flag to know if TLS is being run
  
-;CASE getenv('DUSTEM_WHICH') OF
 CASE !dustem_which OF
   'WEB3p8': BEGIN
     frmt='(A,D,D,D,D,D,I,A)'
@@ -56,7 +96,7 @@ CASE !dustem_which OF
       Nlines=Nlines+1
     ENDWHILE
     close,unit
-free_lun,unit
+    free_lun,unit
     ;== now read the file
     openr,unit,file,/get_lun
     ncurrent=0L
@@ -67,32 +107,21 @@ free_lun,unit
     ENDREP UNTIL first_char NE '#'
     key_str=str
     readf,unit,G0 & ncurrent=ncurrent+1
-    ;===
-;    readf,unit,nspecies
     nspecies=nlines-ncurrent
-;    stop
-;; # grain type, nsize, type keywords, Mdust/MH, rho, amin, amax, alpha/a0 [, at, ac, gamma (ED)] [, au, zeta, eta (CV)]
     one_st={grain_type:'',nsize:0L,type_keywords:'',Mdust_o_MH:0.D0,rho:0.D0,amin:0.D0,amax:0.D0,alpha_o_a0:0.D0, $
             at:0.D0,ac:0.D0,gamma:0.D0, $
             au:0.D0,zeta:0.D0,eta:0.D0}
     st=replicate(one_st,nspecies)
     FOR i=0L,nspecies-1 DO BEGIN
       readf,unit,str
-;VG : remove leading blanks
-      ;strv=str_sep(strcompress(str),' ')
       strv=str_sep(strcompress(strtrim(str,2)),' ')
-;VG
       Nstrv=n_elements(strv)
       ii=0L
       st(i).grain_type=strv(ii) & ii=ii+1
       st(i).nsize=long(strv(ii)) & ii=ii+1
       st(i).type_keywords=strv(ii) & ii=ii+1
-;VG
       IF stregex(st(i).type_keywords, 'pol', /bool) THEN !run_pol = 1
-;VG
-;DP
       IF stregex(st(i).type_keywords, 'dtls', /bool) THEN !run_tls = 1
-
       st(i).Mdust_o_MH=double(strv(ii)) & ii=ii+1
       st(i).rho=double(strv(ii))  & ii=ii+1
       st(i).amin=double(strv(ii)) & ii=ii+1
@@ -117,7 +146,6 @@ free_lun,unit
         st(i).eta=double(strv(ii)) & ii=ii+1
       ENDIF
     ENDFOR
-;    full_st={keywords:key_str,G0:G0,Ngrains:fix(nspecies),grains:st}
     full_st=st
     close,unit
     free_lun,unit
@@ -147,7 +175,6 @@ free_lun,unit
       st(i).ndiscr=fix(strv(6))
       st(i).flag=strv(7)
     ENDFOR
-;    full_st={keywords:key_str,G0:G0,Ngrains:fix(nspecies),grains:st}
     full_st=st
     close,unit
     free_lun,unit
@@ -169,8 +196,9 @@ free_lun,unit
   END
 ENDCASE
  
-RETURN,full_st
+the_end:
  
+RETURN,full_st
  
 END
  
@@ -6,10 +6,15 @@ FUNCTION dustem_read_mix,file,silent=silent
   openr,unit,file,/get_lun
  
 ;==read comments
-  str='' & first_char='#'
+  str=''
+  first_char='#'
+  comments=['']
+  Ncomments=0L
   WHILE first_char EQ '#' DO BEGIN
      readf,unit,str
      first_char=strmid(str,0,1)
+     comments=[comments,str]
+     IF first_char EQ '#' THEN Ncomments=Ncomments+1
   ENDWHILE
  
 ;==read FMIX
@@ -22,7 +27,7 @@ FUNCTION dustem_read_mix,file,silent=silent
   fmix = tmp
  
 ;==save and return structure
-  full_st={file:file,fmix:fmix}
+  full_st={file:file,fmix:fmix,comments:comments}
  
   close,unit
   free_lun,unit
@@ -57,6 +57,8 @@ ENDCASE
 ;spawn,'dustem'
 ;st=dustem_read_all_res(getenv('DUSTEM_RES'),/silent)
 ;print,'params',params
+;stop
+
 spawn,!dustem_f90_exec
  
 st=dustem_read_all_res(!dustem_res,/silent)
 PRO dustem_set_params,params
  
 ;if empty description, do nothing
-  ind=where(*(*!dustem_fit).param_descs EQ '',count)
-
-IF count EQ n_elements(*(*!dustem_fit).param_descs) THEN goto,the_end
+ind=where(*(*!dustem_fit).param_descs EQ '',count)
+IF count EQ n_elements(*(*!dustem_fit).param_descs) THEN BEGIN
+  message,'param_descs not set. Doing nothing',/continue
+  stop
+  goto,the_end
+ENDIF
  
 Nparams=n_elements(*(*!dustem_fit).param_descs)
  
@@ -14,7 +17,7 @@ FOR i=0L,Nparams-1 DO BEGIN
    IF (*(*!dustem_fit).param_func)(i) EQ 0 THEN BEGIN
      str=(*(*!dustem_fit).param_descs)(i)+'=params(i)'
      toto=execute(str)
-     print,params(i)
+     ;print,params[i]
      IF !dustem_verbose NE 0 THEN message,str,/info
    ENDIF ELSE BEGIN
      length=strlen((*(*!dustem_fit).param_descs)(i))
@@ -22,12 +25,9 @@ FOR i=0L,Nparams-1 DO BEGIN
 ;     stop
      IF strr NE 'continuum' THEN BEGIN
        IF strr eq 'isrf2' then strr = 'isrf'
-;       IF getenv('DUSTEM_WHICH') EQ 'VERSTRAETE' THEN BEGIN
        IF !dustem_which EQ 'VERSTRAETE' THEN BEGIN
-;         ddir=GETENV('DUSTEM_DAT')+'les_DAT/'
          ddir=!dustem_dat+'les_DAT/'
        ENDIF ELSE BEGIN
-;         ddir=GETENV('DUSTEM_DAT')
          ddir=!dustem_dat
        ENDELSE
 ;      JPB: Caution this line reads ISRF.DAT through an execute statement
@@ -40,9 +40,7 @@ FOR i=0L,Nparams-1 DO BEGIN
      ENDIF
   ENDELSE
 ENDFOR
-;print,'before dustem_write'
-;stop
-;dustem_write_all,*!dustem_params,getenv('DUSTEM_DAT')
+
 dustem_write_all,*!dustem_params,!dustem_dat
  
 the_end:
@@ -14,38 +14,39 @@ dustem_write_lambda,file_out,st.lambda
  
 ;== GRAIN
 file_out=dir_out_dat+'GRAIN.DAT'
-;file_out=dir_out_dat+((*!dustem_inputs).grain)
 dustem_write_grain_web3p8,file_out,st.grains
  
+;== GAS
+file_out=dir_out_dat+'GAS.DAT'
+dustem_write_gas,file_out,st.gas
+
+stop
+
+;== SPIN
+dustem_write_spin,dir_out_dat,st.spin
+
+;== CHRG
+dustem_write_chrg,dir_out_dat,st.chrg
+
 ;== MIX
-file_out=dir_out_dat
-dustem_write_mix,file_out,st.mix
+dustem_write_mix,dir_out_dat,st.mix
  
 ;== SIZE
-file_out=dir_out_dat
-dustem_write_size_lv,file_out,st.size
+dustem_write_size_lv,dir_out_dat,st.size
  
 ;== QABS
-file_out=dir_out_qabs
-dustem_write_qabs_lv,file_out,st.qabs
+dustem_write_qabs_lv,dir_out_qabs,st.qabs
  
 ;== CALOR
-file_out=dir_out_capa
-dustem_write_calor_lv,file_out,st.calor
+dustem_write_calor_lv,dir_out_capa,st.calor
  
-;VGb : add polarization
 ;== POL
-file_out=dir_out_dat
-dustem_write_pol,file_out,st.pol
-file_out=dir_out_qabs
-for i_axis = 1, 3 do dustem_write_qabspol,file_out,st,i_axis
-;VGe
+dustem_write_pol,dir_out_dat,st.pol
+FOR i_axis = 1, 3 DO dustem_write_qabspol,dir_out_qabs,st,i_axis
  
-;DP : add TLS
+;== TLS
 file_out=dir_out_dat
-;stop
-dustem_write_tls,file_out,st.tls
-file_out=dir_out_qabs
+dustem_write_tls,dir_out_dat,st.tls
  
 sortie:
  
 PRO dustem_write_grain_web3p8,file,st
  
+;Remark: It is a bit odfd that this is using both st and (*!dustem_params)
+;Remark: This should use automated comment detection, as in dustem_write_gas.pro
+
 Ncomments=100
 c=strarr(Ncomments)
  
 PRO dustem_write_mix,dir,st
  
-;Writes ALL the MIX_xxx.DAT files
+IF ptr_valid(st[0]) THEN BEGIN
+  comments=(*st[0]).comments
+  Ncomments=n_elements(comments)
+ENDIF ELSE BEGIN
+  Ncomments=0L
+ENDELSE
  
-  Ncomments=5
-  c=strarr(Ncomments)
+Nst=n_elements(st)
  
-  c(0)='# DUSTEM: mixing coefficient'
-  c(1)='# fraction of neutral PAHs'
-  c(2)='#'
-  c(3)='# f_mix'
-  c(4)='#'
-  
-  Nst=n_elements(st)
-  
-;stop
-  
-  frmt='(E14.6)'
+
+frmts=[ $
+         '(A)', $
+         '(A)', $
+         '(A)', $
+         '(A)', $
+         '(A)']
  
 ;Start loop over type of grain
-  FOR i=0L,Nst-1 DO BEGIN
-;Test if the pointer is pointing to something
-     IF ptr_valid(st(i)) THEN BEGIN
-;Get filename
+FOR i=0L,Nst-1 DO BEGIN
+  ;Test if the pointer is pointing to something
+  IF ptr_valid(st(i)) THEN BEGIN
+        ;Get filename
         ffile=(*st(i)).file
         fv=str_sep(ffile,'/')
         file=dir+fv(n_elements(fv)-1)
-;Open file
+        ;Open file
+        message,'Will write MIX file '+file,/continue
         openw,unit,file,/get_lun
-;Print comments
-        FOR ii=0,Ncomments-1 DO BEGIN
-           printf,unit,c(ii)
-        ENDFOR
-;Print FMIX
+        ;Print comments
+        IF Ncomments NE 0 THEN BEGIN
+          FOR ii=0,Ncomments-1 DO BEGIN
+             printf,unit,comments[ii]
+          ENDFOR
+        ENDIF
+        ;Print FMIX
         FOR k=0L,n_elements((*st(i)).fmix)-1 DO BEGIN
-           printf,unit,((*st(i)).fmix)(k),format=frmt
+            printf,unit,((*st(i)).fmix)(k),format=frmts[k]
         ENDFOR
         close,unit
         free_lun,unit
-     ENDIF
-  ENDFOR
+   ENDIF ELSE BEGIN
+      ;left blank intentionally
+   ENDELSE
+ENDFOR
  
 END