Astronomer: MQL recipes

Recipes

Kepler's third law from planet data
Distance and luminosity of Vega
Luminosity and peak wavelength of Sirius A
Position of Halley's comet with the Kepler equation
Airy disk of a 200 mm telescope
Julian date, sidereal time and Moon age of an observing night
Altitude, azimuth and air mass of Betelgeuse
Sunrise, sunset and equation of time at the winter solstice
Shadow and Hawking temperature of Sagittarius A*
Luminosity distance and lookback time in Lambda-CDM

executed on MentDB

Kepler's third law from planet data

Fit the periods of six planets against their semi-major axes with a power law, check that the exponent is 3/2, and predict the period of Ceres.

Pipeline

01 · dataa in au, T in yearsstat fit→02 · readExponent and r2json select→03 · predictPeriod of Ceresmath eval→04 · convertYears to daysunit convert

Script

admin

#1. Semi-major axis in au and period in years of six planets;
-> "[a]" "[0.387, 0.723, 1.000, 1.524, 5.203, 9.537]";
-> "[T]" "[0.2408, 0.6152, 1.0000, 1.8808, 11.862, 29.457]";
#2. Power law fit T = k*a^b, Kepler expects b = 1.5;
-> "[fit]" (stat fit [a] [T] "power");
json load "rx_kepler" [fit];
-> "[b]" (math round (json select "rx_kepler" /b) 4);
#3. Predict the period of the dwarf planet Ceres at 2.767 au;
-> "[Tc]" (math eval (json select "rx_kepler" /expression) "{\"x\": 2.767}");
concat "Exponent b = " [b] ", r2 = " (math round (json select "rx_kepler" /r2) 8) ", period of Ceres = " (math round [Tc] 3) " yr = " (math round (unit convert [Tc] "yr" "d") 0) " days"

Result

mentdb

Exponent b = 1.4998, r2 = 0.99999977, period of Ceres = 4.603 yr = 1681 days

Go further

Fit the moons of Jupiter the same way and get the mass of Jupiter from k = 4 pi^2 / GM
Draw the log-log fit with plot fit

executed on MentDB

Distance and luminosity of Vega

From the Hipparcos parallax and the apparent magnitude of Vega, get its distance in parsecs and light years, its absolute magnitude through the distance modulus and its luminosity in solar units.

Pipeline

01 · distanced = 1/p in pcmath calc→02 · convertpc to light yearsunit convert→03 · magnitudeDistance modulus and L/L_sunmath calc

Script

admin

#1. Parallax of Vega in milliarcseconds and its apparent magnitude;
-> "[p]" 130.23;
-> "[m]" 0.03;
#2. Distance in parsecs, then in light years;
-> "[d]" (math calc "1000/[p]");
-> "[ly]" (unit convert [d] "pc" "ly");
#3. Distance modulus gives the absolute magnitude, then the luminosity in suns with M_sun = 4.83;
-> "[M]" (math calc "[m]-5*log10([d]/10)");
-> "[L]" (math calc "10^(0.4*(4.83-[M]))");
concat "d = " (math round [d] 3) " pc = " (math round [ly] 2) " ly, M = " (math round [M] 3) ", L = " (math round [L] 1) " L_sun"

Result

mentdb

d = 7.679 pc = 25.04 ly, M = 0.604, L = 49 L_sun

Go further

Propagate the parallax error on the distance with num uncertainty
Add the bolometric correction for hot stars

executed on MentDB

Luminosity and peak wavelength of Sirius A

Use the Stefan-Boltzmann law with the radius and effective temperature of Sirius A, then Wien's law, to get its luminosity and the wavelength of maximum emission.

Pipeline

01 · radiusSolar radii to metresunit calc→02 · luminosity4 pi R^2 sigma T^4math calc→03 · wienPeak wavelength b/T in nmunit calc→04 · roundReadable figuresmath round

Script

admin

#1. Radius and effective temperature of Sirius A;
-> "[R]" (unit calc "1.711*695700 km" "m");
-> "[Teff]" 9940;
#2. Stefan-Boltzmann law L = 4 pi R^2 sigma T^4, compared to the Sun 3.828e26 W;
-> "[L]" (math calc "4*pi*[R]^2*_sigma*[Teff]^4");
#3. Wien law gives the peak wavelength;
-> "[peak]" (unit calc (concat "_b/(" [Teff] " K)") "nm");
concat "L = " (math round (math calc "[L]/1e27") 3) "e27 W = " (math round (math calc "[L]/3.828e26") 1) " L_sun, peak at " (math round [peak] 1) " nm"

Result

mentdb

L = 9.856e27 W = 25.7 L_sun, peak at 291.5 nm

Go further

Get the radius of a star from L and T instead, by inverting the law
Integrate the Planck law with num integrate to get the visible fraction

executed on MentDB

Position of Halley's comet with the Kepler equation

Ten years after perihelion, solve the Kepler equation E - e sin E = M for Halley's comet, then get its true anomaly and its distance from the Sun.

Pipeline

01 · anomalyMean anomaly 2 pi t/Pmath calc→02 · solveKepler equationnum solve→03 · positionTrue anomaly and rmath calc→04 · convertau to kmunit convert

Script

admin

#1. Halley's comet: eccentricity, semi-major axis in au, period in years;
-> "[e]" 0.96714;
-> "[a]" 17.834;
-> "[P]" 75.32;
#2. Mean anomaly 10 years after perihelion;
-> "[M]" (math calc "2*pi*10/[P]");
#3. Solve the Kepler equation E - e sin E = M;
-> "[E]" (num solve (concat "x-" [e] "*sin(x)-" [M]) "x" 0 3.14159);
#4. True anomaly and distance from the Sun;
-> "[nu]" (math calc "2*atan(sqrt((1+[e])/(1-[e]))*tan([E]/2))");
-> "[r]" (math calc "[a]*(1-[e]*cos([E]))");
concat "E = " (math round [E] 5) " rad, true anomaly = " (math round (math rad_to_deg [nu]) 2) " deg, r = " (math round [r] 3) " au = " (math round (math calc (concat (unit convert [r] "au" "km") "/1e6")) 0) " million km"

Result

mentdb

E = 1.78021 rad, true anomaly = 168.05 deg, r = 21.42 au = 3204 million km

Go further

Loop over the dates to build an ephemeris table
Integrate the two-body problem with num ode and compare

executed on MentDB

Airy disk of a 200 mm telescope

The diffraction limit of a telescope comes from the first zero of the Bessel function J1. The script gets the famous factor 1.22, the angular radius of the Airy disk in arcseconds, the energy inside the first dark ring, and the profile near the centre with a Taylor series.

Pipeline

01 · zeroFirst zero of J1num solve→02 · radiusAngular radius in arcsecunit convert→03 · energyEncircled energybessel_j(nu, x)→04 · seriesProfile near the centrecas series

Script

admin

#1. Telescope of aperture D = 200 mm in green light, 550 nm;
-> "[D]" 0.2; -> "[lambda]" 550e-9;
#2. First dark ring of the Airy pattern: first zero of J1;
-> "[x1]" (num solve "bessel_j(1, x)" "x" 3 4.5);
#3. Angular radius theta = x1/pi * lambda/D, in arcseconds;
-> "[theta]" (math calc "[x1]/pi*[lambda]/[D]");
-> "[arcsec]" (unit convert [theta] "rad" "arcsec");
#4. Energy inside the first dark ring: 1 - J0^2 - J1^2;
-> "[ee]" (math calc "1-bessel_j(0, [x1])^2-bessel_j(1, [x1])^2");
#5. Profile near the centre with the CAS;
-> "[near]" (cas series "(2*bessel_j(1,x)/x)^2" "x" 0 4);
concat "{\"x1\": " [x1] ", \"factor_x1_over_pi\": " (math round (math calc "[x1]/pi") 4) ", \"radius_arcsec\": " (math round [arcsec] 3) ", \"encircled_energy\": " (math round [ee] 4) ", \"I_near_centre\": \"" [near] "\"}"

Result

mentdb

{"x1": 3.831705970207512, "factor_x1_over_pi": 1.2197, "radius_arcsec": 0.692, "encircled_energy": 0.8378, "I_near_centre": "1-x^2/4+5/192*x^4"}

Go further

Compare with the seeing: a 1 arcsec seeing hides the diffraction limit of apertures above about 14 cm.
Draw the Airy pattern with plot function "(2*bessel_j(1,x)/x)^2" "x" -10 10.

executed on MentDB

Julian date, sidereal time and Moon age of an observing night

Before pointing the telescope of the Pic du Midi, the observer needs the Julian date of the UTC instant, the local sidereal time (it gives which right ascensions are on the meridian) and the age and illuminated fraction of the Moon, which brightens the sky.

Pipeline

01 · epochUTC instant to Unix epochdate to_epoch→02 · jdJulian date and GMSTmath calc→03 · wrapKeep the fraction of the day and of the lunationmath frac_part→04 · formatHours and minutesmath trunc

Script

admin

#1. Observing night at the Pic du Midi, longitude 0.1423 E, instant given in UTC;
-> "[utc]" "2026-11-01 21:30:00";
-> "[lon]" 0.1423;
#2. Unix epoch of the UTC instant, then Julian date;
-> "[epoch]" (date to_epoch (date from_utc [utc]));
-> "[jd]" (math calc "[epoch]/86400+2440587.5");
#3. Greenwich mean sidereal time, then local sidereal time, in hours;
-> "[gmst]" (math calc (concat "24*" (math frac_part (math calc "(18.697374558+24.06570982441908*([jd]-2451545))/24"))));
-> "[lst]" (math calc (concat "24*" (math frac_part (math calc "([gmst]+[lon]/15)/24"))));
#4. Age of the Moon since the new moon of 2000-01-06 and illuminated fraction;
-> "[age]" (math calc (concat "29.530588853*" (math frac_part (math calc "([jd]-2451550.1)/29.530588853"))));
-> "[illum]" (math calc "(1-cos(2*pi*[age]/29.530588853))/2");
concat "JD = " (math round [jd] 5) ", LST = " (math trunc [lst]) " h " (math round (math calc (concat "60*" (math frac_part [lst]))) 1) " min, Moon age = " (math round [age] 2) " d, illuminated = " (math round (math calc "100*[illum]") 1) " %"

Result

mentdb

JD = 2461346.39583, LST = 0 h 15.3 min, Moon age = 21.67 d, illuminated = 55.1 %

Go further

Use date convert_tz to give the instant in the local time of the observatory.
Replace the mean lunation by the true phase from the elongation of the Moon to check the error of about one day.

executed on MentDB

Altitude, azimuth and air mass of Betelgeuse

From the equatorial coordinates of Betelgeuse and the local sidereal time, the script gets the hour angle, the altitude and the azimuth (quadrant-safe with atan2), then the air mass with the secant of the zenith distance and the atmospheric extinction in magnitudes.

Pipeline

01 · hourHour angle H = LST - RAmath calc→02 · azimuthAzimuth from northmath atan2→03 · airmassAir mass sec zmath sec→04 · degreesRadians to degreesmath rad_to_deg

Script

admin

#1. Betelgeuse J2000: RA 5h55m10.3s, Dec +7 deg 24 min 25 s, seen from Paris;
-> "[ra]" (math calc "15*(5+55/60+10.3/3600)");
-> "[dec]" (math calc "(7+24/60+25/3600)*pi/180");
-> "[lat]" (math calc "48.8566*pi/180");
#2. Local sidereal time read on the observatory clock: 3 h 40 min, hour angle in radians;
-> "[H]" (math calc "(15*(3+40/60)-[ra])*pi/180");
#3. Altitude from the spherical triangle, azimuth from north with atan2;
-> "[alt]" (math calc "asin(sin([lat])*sin([dec])+cos([lat])*cos([dec])*cos([H]))");
-> "[az]" (math calc (concat (math rad_to_deg (math atan2 (math calc "sin([H])") (math calc "cos([H])*sin([lat])-tan([dec])*cos([lat])"))) "+180"));
#4. Air mass with the secant of the zenith distance and extinction for k = 0.20 mag;
-> "[X]" (math sec (math calc "pi/2-[alt]"));
concat "altitude = " (math round (math rad_to_deg [alt]) 2) " deg, azimuth = " (math round [az] 2) " deg, air mass = " (math round [X] 3) ", extinction = " (math round (math calc "0.20*[X]") 3) " mag"

Result

mentdb

altitude = 39.74 deg, azimuth = 134.17 deg, air mass = 1.564, extinction = 0.313 mag

Go further

Use the Kasten-Young formula for the air mass near the horizon, where sec z diverges.
Loop on the sidereal time to get the transit and the best observing window.

executed on MentDB

Sunrise, sunset and equation of time at the winter solstice

For the Observatory of Haute-Provence on 21 December, the script gets the day of the year, the solar declination, the equation of time, the sunset hour angle with refraction and the local times of sunrise and sunset.

Pipeline

01 · dayDay of the yeardate day_of_year→02 · sunDeclination and equation of timemath calc→03 · angletan phi tan deltamath tan→04 · clockHours and minutesmath frac_part→05 · nameDay namedate dayname

Script

admin

#1. Observatory of Haute-Provence, latitude 43.93 N, longitude 5.71 E, time zone UTC+1;
-> "[lat]" (math deg_to_rad 43.93);
-> "[day]" "2026-12-21";
-> "[N]" (date day_of_year [day]);
#2. Solar declination and equation of time in minutes for this day of the year;
-> "[dec]" (math calc "asin(sin(-23.44*pi/180)*cos(2*pi/365*([N]+10)))");
-> "[B]" (math calc "2*pi*([N]-81)/365");
-> "[eot]" (math calc "9.87*sin(2*[B])-7.53*cos([B])-1.5*sin([B])");
#3. Hour angle of sunset with cos w0 = -tan phi tan delta, -0.833 deg for refraction and solar radius;
-> "[w0]" (math calc (concat "acos((sin(-0.833*pi/180)-sin([lat])*sin([dec]))/(cos([lat])*cos([dec])))"));
-> "[noon]" (math calc "12+1-5.71/15-[eot]/60");
#4. Day length, local sunrise and sunset;
-> "[half]" (math calc "[w0]*180/pi/15");
-> "[rise]" (math calc "[noon]-[half]"); -> "[set]" (math calc "[noon]+[half]");
concat (date dayname [day] "en") " " [day] " day " [N] ": declination = " (math round (math rad_to_deg [dec]) 2) " deg, equation of time = " (math round [eot] 2) " min, tan phi tan delta = " (math round (math calc (concat (math tan [lat]) "*" (math tan [dec]))) 4) ", sunrise " (math trunc [rise]) " h " (math round (math calc (concat "60*" (math frac_part [rise]))) 0) " min, sunset " (math trunc [set]) " h " (math round (math calc (concat "60*" (math frac_part [set]))) 0) " min, day length " (math round (math calc "2*[half]") 2) " h"

Result

mentdb

Monday 2026-12-21 day 355: declination = -23.44 deg, equation of time = 1.38 min, tan phi tan delta = -0.4177, sunrise 8 h 9 min, sunset 17 h 3 min, day length 8.89 h

Go further

Add -6, -12 and -18 deg instead of -0.833 to get the civil, nautical and astronomical twilights.
Loop over the 365 days and draw the analemma with plot xy.

executed on MentDB

Shadow and Hawking temperature of Sagittarius A*

The constants of the physics library come as text with their unit, and with their uncertainty when they are measured like G; the script extracts the values with a regular expression, then computes the Schwarzschild radius of the black hole at the centre of the Galaxy, the angular size of its shadow seen by the Event Horizon Telescope and its Hawking temperature.

Pipeline

01 · constantsG, c, hbar, kconstant physics_newtonian_g→02 · extractNumber out of the textstring regexp_extract→03 · distance8.18 kpc in metresunit convert→04 · physicsRs, shadow, Hawking Tmath calc

Script

admin

#1. Constants of the physics library, value extracted from the text;
-> "[re]" "[0-9.]+(E-?[0-9]+)?";
-> "[G]" (string regexp_extract (constant physics_newtonian_g) [re] 0);
-> "[c]" (string regexp_extract (constant physics_c) [re] 0);
-> "[hbar]" (string regexp_extract (constant physics_planck_over_two_pi) [re] 0);
-> "[k]" (string regexp_extract (constant physics_k) [re] 0);
#2. Sagittarius A star: 4.3 million solar masses at 8.18 kpc;
-> "[M]" (math calc "4.3e6*1.98892e30");
-> "[D]" (unit convert 8180 "pc" "m");
#3. Schwarzschild radius and angular diameter of the shadow, sqrt 27 GM/c^2 radius;
-> "[rs]" (math calc "2*[G]*[M]/[c]^2");
-> "[shadow]" (unit convert (math calc "sqrt(27)*[rs]/[D]") "rad" "arcsec");
#4. Hawking temperature hbar c^3 / 8 pi G M k;
-> "[TH]" (math calc "[hbar]*[c]^3/(8*pi*[G]*[M]*[k])");
concat "G = " [G] ", Rs = " (math round (math calc "[rs]/1e9") 2) " million km = " (math round (unit convert [rs] "m" "au") 4) " au, shadow diameter = " (math round (math calc "[shadow]*1e6") 1) " micro arcsec, Hawking T = " (math round_sig (math calc "[TH]*1e14") 3) " x 10^-14 K"

Result

mentdb

G = 6.6743E-11, Rs = 12.7 million km = 0.0849 au, shadow diameter = 53.9 micro arcsec, Hawking T = 1.43 x 10^-14 K

Go further

const info "G" gives the relative uncertainty of G, the least precise fundamental constant.
Do the same for M87*: 6.5 billion solar masses at 16.8 Mpc, the EHT measured 42 micro arcsec.

executed on MentDB

Luminosity distance and lookback time in Lambda-CDM

For a supernova at redshift 1.5 in the Planck 2018 cosmology, the script integrates 1/E(z) to get the comoving and luminosity distances and the distance modulus, then uses the closed form with the inverse hyperbolic sine for the age of the Universe now and at the emission.

Pipeline

01 · hubbleHubble distance and timeunit convert→02 · integralIntegral of dz/E znum integrate→03 · ageAge with asinhmath asinh→04 · modulusDistance modulusmath calc

Script

admin

#1. Flat Lambda-CDM cosmology Planck 2018: H0, matter and dark energy densities, source at z = 1.5;
-> "[H0]" 67.7; -> "[Om]" 0.31; -> "[OL]" 0.69; -> "[z]" 1.5;
#2. Hubble distance c/H0 in Mpc and Hubble time in Gyr;
-> "[DH]" (math calc "299792.458/[H0]");
-> "[tH]" (math calc (concat (unit convert (math calc (concat (unit convert 1 "pc" "km") "*1e6/[H0]")) "s" "yr") "/1e9"));
#3. Comoving distance: integral of dz over E z, then luminosity distance and distance modulus;
-> "[I]" (num integrate (concat "1/sqrt(" [Om] "*(1+x)^3+" [OL] ")") "x" 0 [z]);
-> "[DL]" (math calc "(1+[z])*[DH]*[I]");
-> "[mu]" (math calc "5*log10([DL]*1e6/10)");
#4. Age closed form t = 2/3 tH asinh of sqrt OL/Om times 1+z ^-1.5 / sqrt OL;
-> "[t0]" (math calc (concat "2/3*[tH]/sqrt([OL])*" (math asinh (math calc "sqrt([OL]/[Om])"))));
-> "[tz]" (math calc (concat "2/3*[tH]/sqrt([OL])*" (math asinh (math calc "sqrt([OL]/[Om])*(1+[z])^-1.5"))));
concat "D_C = " (math round (math calc "[DH]*[I]") 0) " Mpc, D_L = " (math round [DL] 0) " Mpc, mu = " (math round [mu] 2) " mag, age now = " (math round [t0] 2) " Gyr, age at z = " (math round [tz] 2) " Gyr, lookback time = " (math round (math calc "[t0]-[tz]") 2) " Gyr"

Result

mentdb

D_C = 4479 Mpc, D_L = 11197 Mpc, mu = 45.25 mag, age now = 13.8 Gyr, age at z = 4.28 Gyr, lookback time = 9.52 Gyr

Go further

Plot D_L against z with plot function and compare with an empty Universe.
Solve with num solve the redshift at which the expansion started to accelerate: Om (1+z)^3 = 2 OL.