OV-GSO-DC / dustem-wrapper_idl

Commit 93b2bb2f4b0b7e25fd036d228a7ec7c42d06dc8f

Authored by Ilyes Choubani
1 parent 690f5f48
Exists in master

other freefree plugin

Showing 1 changed file with 137 additions and 0 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_plugin_ff_hii.pro 0 → 100755
  Diff comments   View file @93b2bb2
... ... @@ -0,0 +1,137 @@
  1 +FUNCTION dustem_plugin_ff_hii,key=key,val=val,scope=scope,paramtag=paramtag,paramdefault=paramdefault,help=help
  2 +
  3 +;+
  4 +; NAME:
  5 +; dustem_plugin_ff_hii
  6 +; PURPOSE:
  7 +; DustEMWrap plugin to compute free-free emission
  8 +; CATEGORY:
  9 +; DustEM, Distributed, Mid-Level, Plugin
  10 +; CALLING SEQUENCE:
  11 +; freefree=dustem_plugin_freefree([,key=][,val=])
  12 +; INPUTS:
  13 +; None
  14 +; OPTIONAL INPUT PARAMETERS:
  15 +; key = input parameter number
  16 +; val = input parameter value
  17 +; OUTPUTS:
  18 +; freefree = free-free spectrum (on dustem wavelengths)
  19 +; OPTIONAL OUTPUT PARAMETERS:
  20 +; scope = scope of the plugin
  21 +; paramdefault = default values of parameters
  22 +; paramtag = plugin parameter names as strings
  23 +; ACCEPTED KEY-WORDS:
  24 +; help = if set, print this help
  25 +; COMMON BLOCKS:
  26 +; None
  27 +; SIDE EFFECTS:
  28 +; None
  29 +; RESTRICTIONS:
  30 +; The DustEMWrap IDL code must be installed
  31 +; PROCEDURE:
  32 +; This is a DustEMWrap plugin
  33 +; EXAMPLES
  34 +; dustem_init
  35 +; vec=dustem_plugin_freefree(scope=scope)
  36 +; MODIFICATION HISTORY:
  37 +; Written by JPB 2022
  38 +; Evolution details on the DustEMWrap gitlab.
  39 +; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
  40 +;-
  41 +
  42 +IF keyword_set(help) THEN BEGIN
  43 + doc_library,'dustem_plugin_ff_hi'
  44 + output=0.
  45 + goto,the_end
  46 +ENDIF
  47 +
  48 +;default values of input parameters
  49 +Te=1.0E4 ;default electron temperature (K)
  50 +;Adding an amplitude just for the sake of it (still not done)
  51 +;Amplitude= ;Amplitude for Te=1.E4 and nu= GHz (arb) and lambda=1E4 microns
  52 +;Not so sure about this emission measure
  53 +EM = 100. ;default emission measure (cm^-6 pc)
  54 +smallp=0.0 ;default polarization fraction
  55 +psi=0. ;default polarization angle
  56 +scope='ADD_SED'
  57 +paramtag=[textoidl('T_{e}')+' [K]','EM','p','Psi [deg]']
  58 +paramdefault=[Te,EM,smallp,psi]
  59 +IF keyword_set(key) THEN BEGIN
  60 + a=where(key EQ 1,count1)
  61 + b=where(key EQ 2,count2)
  62 + c=where(key EQ 3,count3) ;default polarization fraction -newly added
  63 + d=where(key EQ 4,count4) ;default polarization angle -newly added
  64 + IF count1 NE 0 then Te=(val(a))(0)
  65 + IF count2 NE 0 then EM=(val(b))(0)
  66 + IF count3 NE 0 then smallp=val[c[0]]
  67 + IF count4 NE 0 then psi=val[d[0]]
  68 +ENDIF
  69 +
  70 +lambir=dustem_get_wavelengths()
  71 +Nwavs=n_elements(lambir)
  72 +cmic=3.e14
  73 +nu=cmic/lambir ;Hz
  74 +k_B = 1.380649E-23 ;J.K^-1
  75 +g_ff = alog(4.955*1D-2/nu*1D9)+1.5*alog(Te) ;gaunt factor
  76 +tau_ff = 3.014*1D-2*Te^(-1.5)*(nu*1D-9)^(-2)*EM*g_ff ;ff optical depth
  77 +T_ff = Te*(1-exp(-tau_ff))
  78 +
  79 +out=fltarr(Nwavs,3)
  80 +out[*,0] = 2*k_B*T_ff*nu^2/(cmic*1E-6)^2*1E20 ;convert to MJy by multiplying by 1E20
  81 +
  82 +polar_ippsi2iqu,out[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
  83 +
  84 +out[*,1]=Q
  85 +out[*,2]=U
  86 +
  87 +; ;Not so sure about these two:
  88 +; mjy=1 ;output is in MJy/sr
  89 +; lambir_ref=10000.
  90 +;
  91 +;
  92 +; output=fltarr(Nwavs,3) ;newly added
  93 +;
  94 +; use_method='Dickinson2003_norm'
  95 +;
  96 +; CASE use_method OF
  97 +; 'WallsGabaud1998':BEGIN
  98 +; em=1. ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
  99 +; output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
  100 +; fact=Amplitude/I_halpha_R
  101 +; ;=== normalize to the requested Halpha value.
  102 +; output[*,0]=output[*,0]*fact
  103 +; END
  104 +; 'WallsGabaud1998_scaled':BEGIN
  105 +; em=1. ;This is a stupid value for the Emission measure. Result will be rescaled based on I_halpha_R anyway
  106 +; output[*,0]=intensity_free_free(nu,Tgas,em,I_halpha_R=I_halpha_R,nHe=nHe,ergs=ergs,mjy=mjy)
  107 +; norm=interpol(output[*,0],lambir,lambir_ref)
  108 +; output[*,0]=output[*,0]/norm*Amplitude
  109 +; END
  110 +; 'Deschenes2008':BEGIN
  111 +; beta_freefree=2.+1./(10.48+1.5*alog(Tgas/8.e3)-alog(nu/1.e9))
  112 +; output[*,0]=nu^(-1.*beta_freefree)
  113 +; norm=interpol(output[*,0],lambir,lambir_ref)
  114 +; output[*,0]=output[*,0]/norm*Amplitude
  115 +; END
  116 +; 'Dickinson2003_norm':BEGIN
  117 +; ;stop
  118 +; em=1.
  119 +; T4=Tgas/1.e4
  120 +; nu_ghz=nu/1.e9
  121 +; a=0.366*(nu_ghz/1.e9)^(-0.15)*(alog(4.995*1e-2*(nu_ghz)^(-1.))+1.5*alog(Tgas))
  122 +; Tb=8.369e3*a*(nu_ghz)^(-2.)*T4^0.667*10^(0.029*T4)*(1.+0.08) ;in mK for 1 Rayleigh
  123 +; convert_mk_mjy, lambir, Tb, I_Mjy, /RJ
  124 +; norm=interpol(I_Mjy,lambir,lambir_ref)
  125 +; output[*,0]=I_Mjy/norm*Amplitude
  126 +; END
  127 +; ENDCASE
  128 +
  129 +; polar_ippsi2iqu,output[*,0],Q,U,replicate(smallp,Nwavs),replicate(psi,Nwavs)
  130 +;
  131 +; output[*,1]=Q
  132 +; output[*,2]=U
  133 +
  134 +the_end:
  135 +RETURN,out
  136 +
  137 +END
... ...