;+ ; NAME: ; vty16_fig4_2 ; PURPOSE: (one line) ; Plot Young 2016, Fig4-2 ; DESCRIPTION: ; Plot Young 2016, Fig4-2 ; CATEGORY: ; VT3D ; CALLING SEQUENCE: ; vty16_fig4_2 ; INPUTS: ; none ; OPTIONAL OUTPUTS: ; ; SIDE EFFECTS: ; creates vty16_fig_14.png, the plot in the Young 2016. ; vty16_fig3_14_temp_surf.fits - output of vt3s_thermpro ; MODIFICATION HISTORY: ; Written 2013 Jun 5, by Leslie Young, SwRI ; Based on ms_plutoseason/figs/__working_copy/_mimas.pro ; 2016 Mar 24 LAY. Modified for inclusion in vty16 library. ;- pro vty16_fig4_2 ;------------------------- ; this function name ;------------------------- funcname = 'vty16_fig4_2' run = 'vty16_fig4_2' ;======================================= ; SET UP FOR THE CALL ;======================================= ;------------------------- ; Numerical, unit, and physical constants ;------------------------- physconstants constants = { $ ; Numerical constants i : complex(0,1), $ ; sqrt(-1) deg : !dpi/180.,$ ; radian per deg ; Unit conversions hour : 3600.,$ ; s per hour year : 31557600.,$ ; s per year km : 1e5,$ ; cm per km ti_SI : 1.e3, $ ; [cgs thermal inertia] per [SI thermal inertia] ; [J m^-2 K^-1 s^-1/2] per [erg cm^-2 K^-1 s^-1/2] cp_SI : 1.e4, $ ; [cgs specific heat] per [SI specific heat] ; [J kg^-1 K^-1] per [erg g^-1 K^-1] ; Fundemental physical constants, 2002 CODATA values m_u : !phys.m_u, $ ; atomic mass constant, g N_a: 6.0221415d23, $ ; Avogadro constant, mol^-1 k: !phys.k, $ ; Boltzmann constant, erg K^-1 sigma: !phys.sigma, $ ; Stefan-Boltzmann constant erg cm^-2 K^-4 s^-1 G: !phys.g, $ ; gravitational constant, cm^3 g^-1 s^-2 ; Other physical constants sol_norm_1au : 1367.6e3 $ ; solar flux at 1 AU, erg/cm^2/s } ;------------------------- ; Inputs describing the physical problem ;------------------------- ; The physical parameters most often tweaked period = 10. * 24*constants.hour therminertia = 20e3 lat_sol = 0. * constants.deg ; properties of pure N2 dens_n2 = 1.0 ; density of pure alpha N2 at 30 K specheat_n2 = 1.2e7 ; specific heat of pure alpha N2 at 30 K thermcond_n2 = 2.9e4 ; thermal conductivity of alpha N2 at 30 K therminertia_n2 = vt3d_thermalinertia(dens_n2, specheat_n2, thermcond_n2) ; assume therminertia_n2 is split between the density and thermal conductivity dens = dens_n2 * (therminertia/therminertia_n2) specheat = specheat_n2 thermcond = thermcond_n2 * (therminertia/therminertia_n2) inputs = $ { $ ; inputs defining the problem dist_sol_au : 52.45, $ ; ; heliocentric distance in AU albedo : 0.7, $ ; bond albedo, assumed equal to hemispheric emis : 0.90, $ ; emissivity, scalar flux_int : 0.00, $ ; internal heat flux, erg cm^-2 s^-2 therminertia : 5e3, $ ; thermal inertia in erg cm^-2 s^1/2 K^-1, scalar dens : dens, $ ; density g cm^-3, scalar specheat : specheat, $ ; specific heat in erg g^-1 K^-1, scalar thermcond : thermcond, $ ; thermal conductivity in erg cm^-1 s^-1 K^-1, scalar lon_sol_ref : !dpi, $ ; reference sub-solar longitude, radians t_lon_sol_ref : 0.d, $ ; time of lon_sol_ref dlon_sol_dt : -(2.d*!dpi)/period, $ ; rate of change of lon_sol, radian/s lat_sol : lat_sol, $ ; sub-solar latitude is_volatile : 1, $ ; 1 if has a volatile, 0 if bare. flag. scalar mass_volatile : .001, $ ; mass of the volatile slab, g cm^-2, scalar specheat_volatile : 1.3e7, $ ; specific heat of the volatile slab, erg g^-1 K^-1, scalar gravacc : 50., $ ; gravitational acceleration, cm s^-2 name_species : "CH4", $ ; name of the species. String. "N2", "CO", or "CH4" mflux_esc : 0.0 $ ; escape rate in g cm^-2 s^-1 } ;------------------------- ; Grids for the surface, sub-surface, and time ;------------------------- n_per_period = 96L n_period = 5L n_t = n_per_period * n_period zeta = vt3d_zgrid(1.0, /medium) grids = $ { $ lat : 00.0d * constants.deg, $ ; latitude of spot, scalar lon : 0.0d, $ ; longitude of spot, scalar t_start : 0.d, $ ; time of calculation start, s ; t_end : (n_t-1) * period/n_per_period, $ ; time of calculation end, s t_end : n_period * period, $ ; time of calculation end, s tau : 2.d * !dpi/n_per_period, $ ; unitless timestep zeta : zeta $ ; unitless depths of layers } ;------------------------- ; Program control ;------------------------- program = $; : program control { $ n_terms : 7, $ ; number of terms in the sinusoidal expansion iterate_temp_terms : 0, $ ; 1 to iterate the initial temp terms iterate_dfluxdtemp_e : 0, $ ; 1 to iterate the explicit temp terms explicit : 0, $ ; 1 if explicit timestep, 0 if implicit debug : 0 $ ; 1 to print debug statements, 0 otherwise } ;======================================= ; SET UP TO LOOP OVER DISTANCES ;======================================= ;temp_ss = arrgen(30,50,1.) ;dist_sol_au = sqrt(constants.sol_norm_1au * (1-inputs.albedo)/ (inputs.emis * constants.sigma * temp_ss^4)) ; Calculate: ; 2 thermal inertias. ; Typical KBO 20e3 cgs = 20 tiu ; Water 2.1e6 cgs = 2100. tie ; 2 species, N2 and CH4 ; 3 volatile configurations ; Bare ; Local atmosphere ; Local atmosphere and mass slab ; 2 times of day ; dawn, dusk dist_sol_au = arrgen(30.,80, 1.) n_au = n_elements(dist_sol_au) therminertia_arr = [5e3, 2.1e6] & n_therminertia = n_elements(therminertia_arr) ;name_therminertia = ['Low !4C!X', 'High !4C!X'] name_therminertia = ['!4C!X=5 tiu', '!4C!X=2100 tiu'] name_species_arr = ['N2', 'CH4'] n_species = n_elements(name_species_arr) thick_arr = [1, 3] name_is_volatile_arr = ['Bare', 'Local Atm.', 'Local Atm + Slab'] is_volatile_arr = [0, 1, 1] mass_volatile_arr = [0., 0., 1.] vol_linestyle_arr = [0, 2, 3] n_is_volatile = n_elements(name_is_volatile_arr) indx_last = arrgen(n_t-n_per_period, n_t - 1, 1) indx_dawn = indx_last[n_per_period/4] indx_dusk = indx_last[3*n_per_period/4] ;indx_tod = [indx_dawn, indx_dusk] indx_tod = indx_last tod = (indx_tod mod n_per_period)/ float(n_per_period) n_tod = n_elements(indx_tod) tod_color_arr = ['ff0000'xl, '0000ff'xl] temp = fltarr(n_species, n_au, n_therminertia, n_is_volatile, n_tod) p = fltarr(n_species, n_au, n_therminertia, n_is_volatile, n_tod) sol_term_0 = fltarr(n_species, n_au, n_therminertia, n_is_volatile) temp_term_0 = fltarr(n_species, n_au, n_therminertia, n_is_volatile) theta_sva = fltarr(n_species, n_au, n_therminertia, n_is_volatile, 3) freq = abs(inputs.dlon_sol_dt) ; frequency, Hz for i_species = 0, n_species-1 do begin inputs.name_species = name_species_arr[i_species] for i_au = 0, n_au-1 do begin inputs.dist_sol_au = dist_sol_au[i_au] for i_therminertia = 0, n_therminertia-1 do begin inputs.therminertia = therminertia_arr[i_therminertia] for i_is_volatile = 0, n_is_volatile-1 do begin inputs.is_volatile = is_volatile_arr[i_is_volatile] inputs.mass_volatile = mass_volatile_arr[i_is_volatile] vt3d_1loc_diurnal_local, constants, inputs, grids, program, output temp[i_species, i_au, i_therminertia, i_is_volatile, *] = $ output.s[indx_tod].temp_0 ; thermal parameters sol_term_0[i_species, i_au, i_therminertia, i_is_volatile] = output.flux_sol_terms[0] temp_t0 = output.temp_terms[0] temp_term_0[i_species, i_au, i_therminertia, i_is_volatile] = temp_t0 dfluxdtemp_e = vt3d_dfluxdtemp_emit(inputs.emis, temp_t0) dfluxdtemp_s = vt3d_dfluxdtemp_substrate(freq, inputs.therminertia) ; erg cm^-2 s^-1 K^-1 : float dfluxdtemp_v= vt3d_dfluxdtemp_slab(freq, inputs.mass_volatile, inputs.specheat_volatile) dfluxdtemp_a = vt3d_dfluxdtemp_atm(freq, temp_t0, 1.0, inputs.gravacc, inputs.name_species) theta_sva[i_species, i_au, i_therminertia, i_is_volatile,0] = dfluxdtemp_s/(dfluxdtemp_e/4) theta_sva[i_species, i_au, i_therminertia, i_is_volatile,1] = dfluxdtemp_v/(dfluxdtemp_e/4) theta_sva[i_species, i_au, i_therminertia, i_is_volatile,2] = dfluxdtemp_a/(dfluxdtemp_e/4) endfor endfor ; therminertia endfor ; au p[i_species,*,*,*,*] = bz1980_vp(inputs.name_species, temp[i_species,*,*,*,*]) endfor ; species vty16_fig4_2_plot, dist_sol_au, tod, theta_sva, temp_term_0, temp, p, $ name_therminertia, n_per_period end