While I got nostalgic, I looked around on a backup and I found
@hellfire's 2009 Calculate Function that I later incorporated into Lightcone7. For what it's worth, here it is:
function Calculate(OmegaM, OmegaR, OmegaL, HubbleP, zo, Omegat)
{
// >>>>>> Latest update 12-05-2009
// >>>>>> Original:
http://www.geocities.com/alschairn/cc_e.htm
// >>>>>> input
var Ol = OmegaL * 1; // Lambda density par
var Or = OmegaR * 1; // Radiation density par
var Om = OmegaM * 1; // matter density par
var Ok = 1 - Om - Or - Ol; //curvature densioty par
var H = HubbleP * 1; // Hubble const now
var zt = zo * 1; // redshift of object
var at = (1 / (1 + zt)); // requested redshift value
// >>>>>> output
var T = 0; // time
var Dc = 0; // comoving distance
var Dp = 0; // proper distance
var Da = 0; // angular diameter distance
var Dl = 0; // luminosity distance
var Hp = 0; // time variable Hubble constant
var Omt = 0;
var Olt = 0;
var Ort = 0;
// >>>>>> auxiliary
var N = 100000;
var a = 0;
var a_ic = 0;
var qa = 0;
var pa = 1 / N; // Incremental a
var Eset = 0;
var Dtc = 0;
var Dtp = 0;
// >>>>>> conversion 1/H to Myr
var Tyr = 978000;
// >>>>>> Loop from a = 0 to a = 1, stop to get at values
while (a_ic < 1)
{
a = a_ic + (pa / 2); // expansion factor as incremented
qa = Math.sqrt((Om / a) + Ok + (Or / (a * a)) + (Ol * a * a)); // time variable density function (Or input 10000 times hi)
T = T + (pa / qa); // time (age)
Dc = Dc + ((1 / (a * qa)) * pa); // proper distance now
Dp = Dp + ((1 / qa) * pa); // proper distance then
a_ic = a_ic + pa; // new a
if ((a > at) && (Eset == 0))
{
window.document.cc_e.Epoch.value = Math.round((T * (Tyr / H)) * 10000) / 10000;
Dtc = Dc;
Dtp = Dp;
Dtl = Dl;
Dta = Da;
Eset = 1;
};
}
// >>>>>> auxiliary
Okc = Ok;
// >>>>>> Angular diameter distance -- D. Hogg, astro-ph/9905116
if ((Math.abs(Okc)) < 0.05)
{
Da = at * (Dc - Dtc);
}
else if (Okc > 0)
{
Da = at * (1/Math.sqrt(Okc)) * 0.5 * (Math.exp(Math.sqrt(Okc) * (Dc - Dtc)) - Math.exp(- Math.sqrt(Okc) * (Dc - Dtc)));
}
else
{
Okc = - Okc;
Da = at * (1/Math.sqrt(Okc)) * Math.sin(Math.sqrt(Okc) * (Dc - Dtc));
}
if (Da < 0)
{
Da = - Da;
}
// >>>>>> Luminosity distance
Dl = Da / (at * at);
// >>>>>> Conversion
T = T * (Tyr / H);
Dc = Dc * (Tyr / H);
Dtc = Dtc * (Tyr / H);
Dp = Dp * (Tyr / H);
Dtp = Dtp * (Tyr / H);
Da = Da * (Tyr / H);
Dl = Dl * (Tyr / H);
// >>>>>> Output "now"
window.document.cc_e.Age.value = Math.round(T * 10000) / 10000;
window.document.cc_e.Size.value = Math.round(Dc * 10000) / 10000;
window.document.cc_e.Dcomoving.value = Math.round((Dc - Dtc) * 10000) / 10000;
window.document.cc_e.Dproper.value = Math.round((Dp - Dtp) * 10000) / 10000;
window.document.cc_e.Dluminosity.value = Math.round(Dl * 10) / 10;
window.document.cc_e.Dangular.value = Math.round(Da * 10000) / 10000;
//window.document.cc_e.Vrec.value = Math.round((1-at*at) / (1+at*at)*10000)/10000; //temp Vrec calc using normal Doppler
window.document.cc_e.Vrec.value = Math.round(((((Dc - Dtc) / 3.26) * H) / 300000) * 10000) / 10000;
// >>>>>> Hubble parameter then, at a = at
Hp = (1/at) * H * Math.sqrt((Om * (1/at)) + (Or * (1/(at * at))) + (Ok) + (Ol * at * at));
// >>>>>> Omegas
Omt = (Om * (H * H) / (Hp * Hp * at * at * at));
Olt = (Ol * (H * H) / (Hp * Hp));
Ort = (Or * (H * H) / (Hp * Hp * at * at * at * at));
// >>>Ot = (Omt + Olt + Ort);
Ot = 0.0001878754 * Hp * Hp;
//Ot = Ot + ""
//Ot = Ot.substring(0,4)
//parseFloat(Ot)
// >>>>>> auxiliary
Okc = 1 - Omt - Olt - Ort;
// >>>>>> Output "then"
window.document.cc_e.Hthen.value = Math.round(Hp * 10000) / 10000;
window.document.cc_e.Dthen.value = Math.round((Dc - Dtc) * at * 10000) / 10000;
window.document.cc_e.Vrecthen.value = Math.round(((((Dc - Dtc) * at / 3.26) * Hp) / 300000) * 10000) / 10000;
window.document.cc_e.Sizethen.value = Math.round(Dtc * 10000) / 10000;
window.document.cc_e.OmegaMt.value = Math.round(Omt * 10000000) / 10000000;
window.document.cc_e.OmegaLt.value = Math.round(Olt * 10000000) / 10000000;
window.document.cc_e.OmegaRt.value = Math.round(Ort * 10000000) / 10000000;
window.document.cc_e.Omegat.value = Math.round(Ot * 10000000) / 10000000;
// window.document.cc_e.Omegat.value = Ot;
}
