;+ ; NAME: ; vty16_fig3_3 ; PURPOSE: (one line) ; Plot Young 2016, Fig3.3 ; DESCRIPTION: ; Plot Young 2016, Fig3.3 ; CATEGORY: ; VT3D ; CALLING SEQUENCE: ; vty16_fig3_3 ; INPUTS: ; none ; OPTIONAL OUTPUTS: ; phase_num - phase for the numeric calculation, radian ; temp_num - surface temperature from the numeric calculation, K ; phase - phase for the approximations, radian ; temp_1_zeta0 - approximate surface temperature, M=1, K ; temp_7_zeta0 - approximate surface temperature, M=7, K ; temp_7_iter_zeta0 - approximate surface temperature, ; M=7 with iterated T0, K ; SIDE EFFECTS: ; creates vty16_fig3_3.png, the plot in the Young 2016. ; vty16_fig3_3_table1.txt - output of vt3s_thermpro ; vty16_fig3_3_table2.txt - approximations ; vty16_fig3_3_table3.txt - temperature terms ; MODIFICATION HISTORY: ; Written 2012 Sep 29, by Leslie Young, SwRI ; 2016 Mar 24 LAY. Modified for inclusion in vty16 library. ;- pro vty16_fig3_3, phase_num, temp_num, phase, temp_1_zeta0, temp_7_zeta0, temp_7_iter_zeta0 ;------------------------- ; this function name ;------------------------- funcname = 'vty16_fig3_3' ;------------------------- ; set up for correct plotting ;------------------------- loadct, 0 device,/decompose !p.background = 'ffffff'xl !p.color=0 !p.charsize=1.5 ;------------------------- ; constants ;------------------------- physconstants deg = !pi/180. sol_norm_1au = 1367.6e3 ; solar constant at 1 AU, erg/cm^2/s s_per_day = 24. * 3600. s_per_year = 365.25 * s_per_day cm_per_au = 149597870.691e5 ; cm per AU km_per_au = 149597870.691 ; km per AU _i = complex(0,1) ;------------------------- ; formats for printing ;------------------------- formf = '("' + funcname + ': ", A-20,"=",F15.6," ",A-10)' ;------------------------- ; Insolation on an asteroid ;------------------------- period = .9424218 * s_per_day n_phase = 720 ; number per period phase = dindgen(n_phase) * 2.d * !dpi / n_phase lat = 30. * deg ; latitude in radians lon = 0. ; east longitude in radians lat_sol = replicate(2.24*deg, n_phase) ; subsolar latitude in radians as a function of time ; subsolar east longitude in radians; decreases with time for ; prograde objects lon_sol = (!pi/2 - phase) ; Calculate the exact insolation ; mu0 is the cosine of the solar illumination angle mu0 = vt3d_solar_mu(lat, lon, lat_sol, lon_sol) albedo = 0.6 dist_sol_au = 9.487299 flux_sol_ss = (1-albedo) * sol_norm_1au / (dist_sol_au^2) flux_sol = flux_sol_ss * reform(mu0) ; Calculate the sinusoidal expansion h_phase0 = lon - lon_sol[0] ; flux_sol_t_1 and flux_sol_t_7 are the M=1 and M=7 arrays of coefficient flux_sol_t_1 = vt3d_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol[0], 1) flux_sol_t_7 = vt3d_sol_terms_diurnal(dist_sol_au, albedo, lat, h_phase0, lat_sol[0], 7) ; flux_sol_1 and flux_sol_7 are the M=1 and M=7 approximations flux_sol_1 = vt3d_solwave(flux_sol_t_1,phase) flux_sol_7 = vt3d_solwave(flux_sol_t_7,phase) ; Calculate the temperature terms flux_int = 0. emis = 1. freq = 2.d * !dpi / period therminertia = 16e3 ; Mimas region 1 dens = 1.0 specheat = 3.e6 thermcond = therminertia^2/(dens * specheat) z_skin = vt3d_skindepth(dens, specheat, thermcond, freq) temp_eq = (flux_sol / (!phys.sigma * emis) ) ^ 0.25 temp_terms_1 = vt3d_temp_terms_bare(flux_sol_t_1,flux_int,emis,freq,therminertia) temp_terms_1_iter = vt3d_temp_terms_bare_iter(flux_sol_t_1,flux_int,emis,freq,therminertia, thermcond) temp_terms_7 = vt3d_temp_terms_bare(flux_sol_t_7,flux_int,emis,freq,therminertia) temp_terms_7_iter = vt3d_temp_terms_bare_iter(flux_sol_t_7,flux_int,emis,freq,therminertia, thermcond) dtemp_dzeta = 0. temp_1_zeta0 = vt3d_thermwave(temp_terms_1, phase, dtemp_dzeta, z_skin, 0.) temp_7_zeta0 = vt3d_thermwave(temp_terms_7, phase, dtemp_dzeta, z_skin, 0.) temp_7_iter_zeta0 = vt3d_thermwave(temp_terms_7_iter, phase, dtemp_dzeta, z_skin, 0.) temp_0 = vt3d_temp_term0_bare(float(flux_sol_t_1[0]), flux_int, emis) phi_e = vt3d_dfluxdtemp_emit(emis, temp_0) phi_s = vt3d_dfluxdtemp_substrate(freq, therminertia) print, 'Theta', 4 * phi_s / phi_e, form=formf ; Call vt3d_thermpro ; This has the same inputs and outputs and John Spencer's thermpro. vt3d_thermpro,tsurf_vt3d,tod_vt3d,tdeep_vt3d,teq_vt3d,balance_vt3d,tdarr,zarr0, $ rhel=dist_sol_au, $ ; heliocentric distance, AU alb=albedo, $ ; bolometric albedo ti=therminertia, $ ; thermal inertia, erg-cgs. rot=period/s_per_day, $ ; Rotation period, days ; Optional physical keyword parameters: emvty=emis, $ ; emissivity (default 1.0) rho=dens, $ ; density, g cm-3. Can be an array, giving the density of each slab. (default 1.0) cp=specheat, $ ; specific heat, erg g-1 K-1 (default H20) plat=lat, $ ; latitude (radians) (default 0.0) sunlat=lat_sol[0], $ ; Subsolar latitude (radians) (default 0.0) heatflow=flux_int, $ ; $ ; Endogenic heat flow, erg cm-2 s-1 (default 0.0) ; Optional simulation keyword parameters: ntinc=1000, $ ; Time increments per day (default 5000) nday=20, $ ; Number of days per run (default 4) ; Optional program control silent = 1 ; tod = time of day, which is zero at midnight. ; At midnight, the hour angle of the sun in pi. ; Add pi to turn it into hour angle, ; then subtract h_phase0 to turn it into phase reorder_by_phase, tod_vt3d-h_phase0+!pi, tsurf_vt3d, phase_num, temp_num window, 0, xs=840, ys=514 !p.color = 0 plot, phase/deg, temp_eq, /nodata, $ xtit='Phase, deg', xstyle=3, xticks=4,xminor=6,$ ytit='Temperature, K', yrange=[60,92], /ystyle ; thick gray = numeric oplot, phase_num/deg, temp_num, co='a0a0a0'xl, th=6 oplot, 270-[5,45], 90+.4+[0,0], color='a0a0a0'xl, thick=6 xyouts, 270, 90, 'Numerical' ; solid M=1 oplot, phase/deg, temp_1_zeta0, line=0, thick=2 oplot, 270-[5,45], 88+.4+[0,0], line=0, thick=2 xyouts, 270, 88, 'Approximation, M=1' ; dot-dot-dashed M=7 oplot, phase/deg, temp_7_zeta0, line=4, thick=2 oplot, 270-[5,45], 86+.4+[0,0], line=4, thick=2 xyouts, 270, 86, 'Approximation, M=7' ; dash M=7, iter oplot, phase/deg, temp_7_iter_zeta0, line=2, thick=2 oplot, 270-[5,45], 84+.4+[0,0], line=2, thick=2 xyouts, 270, 84, 'Approximation, M=7,' xyouts, 270, 82.5, 'iterated T!D0!N' tv2im, funcname + '.png', /corner forprint, phase_num/deg, temp_num, $ TEXTOUT = funcname + '_table1.txt', $ comment = string('phase', 'numeric', for='(A7, 1A9)'), $ FORMAT = '(F7.2, 1F9.4)' forprint, phase/deg, temp_1_zeta0, temp_7_zeta0, temp_7_iter_zeta0, $ TEXTOUT = funcname + '_table2.txt', $ comment = string('phase', 'M=1', 'M=7', 'M=7,iter', for='(A7, 3A9)'), $ FORMAT = '(F7.2, 3F9.4)' pad = replicate(complex(0,0),6) forprint, indgen(8), [temp_terms_1, pad], [temp_terms_1_iter,pad], temp_terms_7, temp_terms_7_iter,$ TEXTOUT = funcname + '_table3.txt', $ comment = string('m', 'M=1', 'M=1,iter', 'M=7', 'M=7,iter', for='(A1, 4A18)'), $ FORMAT = '(I1, 4(F7.2," +",F7.2," i"))' end