pro DustAndIce, mass, diam, dump COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Input Syntax Verification ----------------------------------------------------------- ;-31 mass={val:5.95e18, unc:0.71e18} diam={val:187., unc:18} ;-107 ; mass={val:1.10e19, unc:0.12e19} ; diam={val:254., unc:18} ;-CG ; mass={val:998.2e9, unc:.3e9} ; diam={val:3.299 /2 , unc:0.05/2} help, mass, diam, dump diaMeter={val:diam.val*1000., unc:diam.unc*1000.} vol = {val:(!PI*diaMeter.val^3)/6., unc:3*diam.unc/diam.val} dens = {val: mass.val/vol.val, unc: 0. } dens.unc = dens.val* sqrt( (mass.unc/mass.val)^3 + (3*diam.unc/diam.val)^3. ) help, dens ;tbd: plus sux pour mass/diam/vol pour la bonne distro d'erreur ;sur vol et pas 3dD/D necessairement ; dens = {val:1280., unc:134.} ; met={val:2600., unc:200.} ;-??? met={val:2250., unc:200.} macro = {val:100. *( 1. - dens.val/met.val ), unc:0.} macro.unc = sqrt(2)*macro.val*sqrt( (met.unc/met.val)^2. + (dens.unc/dens.val)^2. ) print, 'dens' print, dens print, 'macro' print, macro ;---camilla ice = 940 ; kg/m3 Si = {min: 2200., $ med: 2600., $ max: 3000 } poro = {min:10., step:20., N:3} densRange=[1000,3500] fracRange=[0,1] ratioRange=[0.1,100] ;--- daphne ice = 940 ; kg/m3 Si = {min: 2000., $ med: 2200., $ max: 2400 } poro = {min:0., step:10., N:3} densRange=[1000,3500] fracRange=[0,1] ratioRange=[0.1,100] ;-CG 67P ; densRange=[1000,3500] ; fracRange=[0,.4] ; ratioRange=[0.1,20] charPlot = 1.1 charText = 1.0 colValid = 'Dark Goldenrod' thkValid = 10 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Porosity and Dust Fraction --------------------------------------------------------- nbFrac = 100 fracDust = indgen(nbFrac)/(nbFrac-1.) pArr = (poro.min + poro.step*indgen(poro.N) )/100. dustDens= fltarr(nbFrac, poro.N) statDens= fltarr(3,poro.N) dustToIceMass = fltarr(nbFrac, poro.N) dustToIceStat = fltarr(3,poro.N) massDust = fltarr(nbFrac, poro.N) massIce = fltarr(nbFrac, poro.N) volDust = fltarr(nbFrac, poro.N) volIce = fltarr(nbFrac, poro.N) volVoid = fltarr(nbFrac, poro.N) rCore = fltarr(nbFrac, poro.N) statMassDust = fltarr(3, poro.N) statMassIce = fltarr(3, poro.N) statVolDust = fltarr(3, poro.N) statVolIce = fltarr(3, poro.N) statVolVoid = fltarr(3, poro.N) statCore = fltarr(3, poro.N) for k=0, poro.N-1 do begin dustDens[*,k] = (dens.val-ice*(1-fracDust-pArr[k]))/fracDust trash = min( abs(dustDens[*,k]-Si.min), pMin) statDens[0,k] = fracDust[pMin] trash = min( abs(dustDens[*,k]-Si.med), pMed) statDens[1,k] = fracDust[pMed] trash = min( abs(dustDens[*,k]-Si.max), pMax) statDens[2,k] = fracDust[pMax] volDust[*,k] = fracDust*vol.val volIce [*,k] = (1-fracDust-pArr[k])*vol.val volVoid[*,k] = pArr[k]*vol.val rCore[*,k] = (diam.val/2.) * fracDust^(1./3.) massDust[*,k] = volDust[*,k] * dustDens[*,k] massIce [*,k] = volIce [*,k] * ice statVolDust[*,k] = volDust[[pMin,pMed,pMax],k] statVolIce [*,k] = volIce [[pMin,pMed,pMax],k] statVolVoid[*,k] = volVoid[[pMin,pMed,pMax],k] statMassDust[*,k] = massDust[[pMin,pMed,pMax],k] statMassIce [*,k] = massIce [[pMin,pMed,pMax],k] statCore[*,k] = rCore[[pMin,pMed,pMax],k] dustToIceMass[*,k] = massDust[*,k]/massIce[*,k] dustToIceStat[*,k] = dustToIceMass[[pMin,pMed,pMax],k] endfor ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- --------------------------------------------------------------------------------- cgPS_open, filename=dump+'.png', /Metric, /Decomposed, /Portrait, $ xsize=15, ysize=20, Language_Level=2, /Quiet xBot = 0.15 yBot = 0.07 xLen = 0.82 yLen = 0.45 ySep = 0.01 ;--III.2-- rosity and Dust Fraction --------------------------------------------------------- pos = [ xBot, yBot+yLen+ySep, xBot+xLen, yBot+2*yLen+ySep ] cgPlot, 0,0, /NoData, /Normal, $ xStyle=1, xRange=fracRange, $ yStyle=1, yRange=densRange, $ yTitle='Dust density (kg.m!U3!N)', $ xTickName=replicate(' ',20), $ position=pos, $ charSize=charPlot cgColorFill, fracRange[[0,1,1,0,0]], $ [Si.min, Si.min, Si.max, Si.max, Si.min], color='Gray', /Data cgPlot, /OverPlot, fracRange, [1,1]*Si.med, lineStyle=2, color='Dark Gray' cgText, 0.025, Si.min+20, 'Organics' , charSize=charText cgText, 0.025, Si.max+20, 'Silicates', charSize=charText cgText, 0.025, Si.med+20, 'Average' , charSize=charText for k=0, poro.N-1 do begin cgPlot, /OverPlot, fracDust, dustDens[*,k] valid = where( dustDens[*,k] ge Si.min and dustDens[*,k] le Si.max ) if valid[0] ne -1 then begin cgPlot, /OverPlot, fracDust[valid], dustDens[valid,k], thick=thkValid, color=colValid endif cgPlot, /OverPlot, statDens[1,k], Si.med, pSym='Filled Circle' cgPlot, /OverPlot, statDens[0,k], Si.min, pSym='Open Diamond', symSize=1.35 cgPlot, /OverPlot, statDens[2,k], Si.max, pSym='Open Square' xKey = 0.90 trahsh = min( abs(fracDust-xKey), pKey ) cgText, xKey, dustDens[pKey[0],k]+20, $ strTrim(string(100*pArr[k],format='(I)'),2)+'%', charSize=charText endfor cgPlot, 0,0, /NoData, /Normal, /NoErase, $ xStyle=1, xRange=fracRange, $ yStyle=1, yRange=densRange, $ position=pos, $ xTickName=replicate(' ',20), $ charSize=charPlot ;--III.2-- rosity and Dust Fraction --------------------------------------------------------- pos = [ xBot, yBot, xBot+xLen, yBot+yLen ] cgPlot, 0,0, /NoData, /NoErase, /Normal, /yLog, $ xStyle=1, xRange=fracRange, $ yStyle=1, yRange=ratioRange, $ xTitle='Dust fraction (volume)', $ yTitle='Dust-to-ice mass ratio', $ position=pos, $ charSize=charPlot for k=0, poro.N-1 do begin cgPlot, /OverPlot, fracDust, dustToIceMass[*,k] ; cgPlot, /OverPlot, statDens[*,k], dustToIceStat[*,k], psym=2 valid = where( dustDens[*,k] ge Si.min and dustDens[*,k] le Si.max ) if valid[0] ne -1 then begin cgPlot, /OverPlot, fracDust[valid], dustToIceMass[valid,k], thick=thkValid, color=colValid endif cgPlot, /OverPlot, statDens[1,k], dustToIceStat[1,k], pSym='Filled Circle' cgPlot, /OverPlot, statDens[0,k], dustToIceStat[0,k], pSym='Open Diamond', symSize=1.35 cgPlot, /OverPlot, statDens[2,k], dustToIceStat[2,k], pSym='Open Square' xKey = 0.10 trahsh = min( abs(fracDust-xKey), pKey ) cgText, xKey, dustToIceMass[pKey[0],k], $ strTrim(string(100*pArr[k],format='(I)'),2)+'%', $ charSize=charText, align=1 endfor cgPS_close, /png, Delete_PS=0 forprint, pArr*100, $ 100*statDens[0,*], 100*statDens[1,*], 100*statDens[2,*], $ ;-frac dust dustToIceStat[0,*], dustToIceStat[1,*], dustToIceStat[2,*], $ ;mass ice-to-dust dustToIceStat[0,*]*ice/Si.med, $ dustToIceStat[1,*]*ice/Si.med, $ dustToIceStat[2,*]*ice/Si.med, $ ;volume ice-to-dust format='(I3,2x,3(2x,I4),2x,3(2x,F5.1),2x,3(2x,F5.1))', $ /Silent end