%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This README file describes some of the files in the source code directory used for integrating the Williams-Boggs orbital dynamics equations, along with a number of "test" and "check" codes. The code was written by Brian Arbic, based upon equations written by Jim Williams, as referenced in Daher et al. (2021). Arbic frequently communicated with Jim Williams as the code was written, and Dale Boggs also performed some imporant double-checks of the code. UPDATE: This README file is specifically for "v9" codes, the code version that was re-run in June 2021, for the revised version of the Daher et al. paper. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The most important change between earlier "v8" and "v7" codes is the change-over to "tidesinMoon2" codes for the impacts of tides in the Moon on da/dt, de/dt, and di/dt. The "tidesinMoon2" codes contain higher-order terms that are not present in the original "tidesinMoon" codes. We also tested "tidesinMoon3" codes that contain even more terms, but Jim was not as confident in the de/dt equation for "tidesinMoon3"; thus I decided to go with the "tidesinMoon2" codes for "v8". % Many of the source code and Matlab output files have a "_nodeLF" % version that includes/writes out the diurnal and semidiurnal % node and L+F constituents. The "nodeLF" codes were added in % May/June 2020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Code for initializing Schindelegger results call_initialize_Schindelegger_results.m calls on initialize_Schindelegger_results.m to invert dissipations to solve for ksinchi values, for the experiments done at a set of discrete rotation rates puts out initialize_Schindelegger_results.mat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Code for setting parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set_parameters_presentday_v9_use_structures.m set_parameters_presentday_v9_use_structures_nodeLF.m -- Sets parameters for present-day. Output files are "parameters_presentday_v9_use_structures.mat" and "parameters_presentday_v9_use_structures_nodeLF.mat" These files store parameter values (such as G, M_E, etc.), constituent index parameters such as the jM, jW, jOmega integers, present-day values of orbital parameters such as a, eccentricity, etc., and present-day values of ksinchi for the body and ocean tides, in four structures The output .mat files are called on several files below and by the "main" files in the run directory The code interp_ocean_ksinchi_presentday_nodeLF.m was written to test the interpolation of k sin chi values from other constituents onto the node and L+F frequencies--this code is not called by any other codes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Other important codes that are called upon by "higher" codes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute_S_E.m -- computes S_E from n and nprime compute_S_A.m -- computes S_A from n and nprime compute_S_W.m -- computes S_W from omega_E and other parameters compute_C_of_omega_E.m -- computes moment of inertia C as a function of omega_E and other parameters compute_dnda.m -- computes dnda via equation (32) of Williams and Boggs (2016). Used to compute dn/dt from da/dt. Strictly for checking numbers in their Table 7; not used elsewhere in code. compute_n_fromaandnprime.m -- computes n from a and nprime, via equation (31) of Williams and Boggs (2016). Wsed when we timestep, to get updated values of n which are then used in some of the d/dt equations. compute_precession_rate.m -- computes precession rate % Codes for computing the "Ucombinations" in the Williams notes: computeUcombinations_diurnaltides.m computeUcombinations_diurnaltides_nodeLF.m computeUcombinations_diurnaltides_obliquity.m computeUcombinations_diurnaltides_obliquity_nodeLF.m computeUcombinations_longperiodtides.m computeUcombinations_longperiodtides_obliquity.m computeUcombinations_semidiurnaltides.m computeUcombinations_semidiurnaltides_nodeLF.m computeUcombinations_semidiurnaltides_obliquity.m computeUcombinations_semidiurnaltides_obliquity_nodeLF.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Codes for computing d/dt of omega_E, epsilon, a, eccentricity, and inclination% due to tides on Earth % (equations 28-30 and 33-34 of Williams and Boggs 2016). Note that the % computations are usually divided up into tidal species: compute_da_dt_diurnal.m compute_da_dt_longperiod.m compute_da_dt_semidiurnal.m compute_deccentricity_dt_diurnal.m compute_deccentricity_dt_longperiod.m compute_deccentricity_dt_semidiurnal.m compute_depsilon_dt_diurnal.m compute_depsilon_dt_longperiod.m compute_depsilon_dt_semidiurnal.m compute_dinclination_dt_diurnal.m compute_dinclination_dt_longperiod.m compute_dinclination_dt_semidiurnal.m compute_domega_E_dt_diurnal.m compute_domega_E_dt_semidiurnal.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Codes used in computing the effects of the "tidesinMoon" compute_noderate.m compute_xi_Moon_coremantledissipation.m--this code is superceded by compute_KoverC_and_xi_lunarCMB.m, which was created in April 2021 in order to simultaneously compute KoverC and xi, prompted by reviewer comments % Codes added in mid-to-late summer 2020, as I moved towards using the % Ward 1975 and Williams et al. 2001 equations for sin I The test code compute_lunarequatorialtilt_asfunctionof_semimajoraxis.m was written to test the Ward 1975 equation for sin I It calls upon the Williams approximation code for comparison tests, but again the Williams appoximation is not used in the "v8" time-stepping codes. The following test codes (not called upon by the main time-stepping codes) compute sin I from the Ward 1975 and Williams et al. 2001 equations, for fixed values of inclination (from the present-day and from the mean of the Monte Carlo simulation values). In the "v8" we had an alternative solution of the Williams et al. 2001 equation described below, but we have taken that out of the v9 codes test_WardEquation_largerinclination.m--computes I as a function of a, from the Ward 1975 equation, using a fixed value of inclination equal to the mean inclination from the main Monte Carlo simulations produces test_WardEquation_largerinclination.mat test_Williams2001_largerinclination.m--ditto, but uses the Williams et al. 2001 equation produces test_Williams2001_largerinclination.mat test_WardEquation_presentdayinclination.m--computes I as a function of a, from the Ward 1975 equation, using a fixed value of inclination equal to the present-day inclination produces test_WardEquation_presentdayinclination.mat test_Williams2001_presentdayinclination.m--ditto, but uses the Williams et al. 2001 equation produces test_Williams2001_presentdayinclination.mat These two codes are called upon by the test codes above, and are also called upon repeatedly by the main "v9" timestepping codes: compute_lunarequatorialtilt_WardEquation.m--computes I from equation in Ward (1975) compute_lunarequatorialtilt_Williams2001.m--computes I from equation (52c) in Williams et al. (2001); calls on root2d_for_Williams2001_variableKoverC.m--an extra program had to be created in order to solve for I in the Williams et al. 2001 equations under the complexity of a variable (non-constant) KoverC_Moon_coremantledissipation % The next three codes compute da/dt, de/dt, di/dt due to tides within % the Moon (and CMB effects) based upon earlier, simpler equations for % the "tidesinMoon" part. compute_da_dt_tidesinMoon.m compute_deccentricity_dt_tidesinMoon.m compute_dinclination_dt_tidesinMoon.m % The next three codes compute da/dt, de/dt, di/dt due to tides within % the Moon (and CMB efects) based upon the higher-order series that % Jim came up with in November 2020. These "tidesinMoon2" codes are % used in the "v8" and "v9" orbital dynamics codes. compute_da_dt_tidesinMoon2.m compute_deccentricity_dt_tidesinMoon2.m compute_dinclination_dt_tidesinMoon2.m % The next three codes compute da/dt, de/dt, di/dt due to tides within % the Moon (and CMB efects) based upon even higher-order series that % Jim came up with in November 2020. These "tidesinMoon3" codes were % tested but not used in the "v8" or "v9" orbital dynamics codes. If someday % they are used they will have to be updated, as the "tidesinMoon2" % codes were, to include the different values of "Q" for different % tidal frequencies within the Moon. compute_da_dt_tidesinMoon3.m compute_deccentricity_dt_tidesinMoon3.m compute_dinclination_dt_tidesinMoon3.m % The codes that I used to test the "tidesinMoon2" and "tidesinMoon3" % codes are test_Williams_equations_January2021.m-->changed to test_Williams_equations_June2021.m in order to employ the "v9" and variable KoverC updates test_Williams_equations_November2020_rewritearguments.m The first code tests the d/dt rates for "tidesinMoon2" and "tidesinMoon3" vs. the original "tidesinMoon" codes--with the new "v8" parameters of January/February 2021. The second code is written to make sense of the subscripts "F", "l", etc. having to do with tidal frequencies within the Moon. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Code used to check present-day frequencies: check_presentday_tidalconstituentfrequencies.m % Code above calls on several other codes above plus this new code: compute_longitudeofperigeerate.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Code used to check ratios of present-day "U" factors vs. ratios of tidal forcing amplitudes taken from Table 1 of Arbic et al. 2004, or in the cases of the diurnal and semidiurnal node and L+F constituents, the new table that Jim Williams sent in June 2019 check_relationship_Ucombinations_to_tidalamplitudes.m check_relationship_Ucombinations_to_tidalamplitudes_nodeLF.m % Codes call on the "computeUcombinations" codes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Now we get to some of the higher-level codes used to check the numbers in Table 7 of Williams and Boggs (2016) % In 2020 we checked our results against the new Table that Jim Williams sent % in a 2019 email %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % We compute the d/dt rates from the codes above (i.e. % "compute_depsilon_dt_diurnal.m", etc.) The values put out by my % codes are in MKS units, and are then converted into units used in % Williams and Boggs Table 7 and Jim's new table, in order to compare. check_presentday_da_dt.m check_presentday_deccentricity_dt.m check_presentday_depsilon_dt.m check_presentday_dinclination_dt.m check_presentday_domega_E_dt.m % This code computes the dinclination/dt, deccentricity/dt, da/dt % values due to the tides and core\mantle dissipation in the Moon: check_presentday_impactoftidesinMoon.m -- changed slightly on 16 % June 2019, in response to 12 June 2019 email from Jim Williams % Changed again on December 24, 2020, February 8, 2021, and June 12, 2021, % in order to incorporate the "tidesinMoon2" codes in the test, and, % in the latter case, in order to incorporate variable KoverC %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % These codes check constituent dissipation values, given values of % ksinchi, and vice versa (with the latter serving as a test of the % inversion of Jim's power dissipation formulae). The M2 body tide % dissipation is compared to the value in Ray et al. (2001), and the M2, % S2, N2, K2, K1, O1, P1, Q1 ocean tide dissipations are compared to % values in Egbert and Ray (2003). %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check_presentday_diurnal_dissipation_and_ksinchi.m check_presentday_longperiod_dissipation_and_ksinchi.m check_presentday_semidiurnal_dissipation_and_ksinchi.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The dissipation/ksinchi codes above call on the codes below (power dissipation equations % P20, P21, and P22 in Williams 6 November 2017 writeup): %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute_diurnal_powerdissipation_from_ksinchi.m compute_longperiod_powerdissipation_from_ksinchi.m compute_semidiurnal_powerdissipation_from_ksinchi.m compute_diurnal_ksinchi_from_powerdissipation.m compute_longperiod_ksinchi_from_powerdissipation.m compute_semidiurnal_ksinchi_from_powerdissipation.m % Time-stepping the orbital evolution equations using results from our % tidal simulations requires us to invert the Williams equations so % that we can obtain "ksinchi" values from power dissipation values. % We will do this inversion with the 3 % "compute...ksinchi_from_powerdissipation.m" codes above. % In turn we use the 3 "compute...powerdissipation_from_ksinchi.m" codes above % to compute power dissipation from body and ocean tides on Earth %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% These codes compute dissipation due to tides on the Moon, and due to the core-mantle boundary effects on the Moon: compute_powerdissipation_tidesinMoon.m -- rewritten 24 December 2020, to account for the fact that dissipation is simply proportional to da/dt compute_powerdissipation_Moon_coremantleboundary.m These codes are not called upon by the main time-stepping routines. This code computes the present-day power dissipation due to tides in the Moon and core-mantle boundary effects: check_dissipation_tidesinMoon_lunarcoremantleboundary.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The codes below compute all the U factors and arranges them into % vectors, to be used in the main time-stepping routines. This code % calls on all six of the "computeUcombinations_" codes listed above: function_compute_all_Ufactors_fortimestepping.m function_compute_all_Ufactors_fortimestepping_nodeLF.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%