0
$\begingroup$

(* this function yields the length of sunrise/sunset in seconds on the nth day
of the year (0.5 <= n <= 366.5) for latitude lat, where lat is measured
in radians *)

temp1942[n_, lat_] = 

6875.493541569879*
 (2.*ArcCos[-0.014543897651582656*Sec[lat]*Sec[0.005782961777094692 - 
        0.4001419318234436*Cos[0.017167172970436028*n] - 
        0.0060922154967620835*Cos[0.034334345940872056*n] - 
        0.002387468786938206*Cos[0.05150151891130809*n] + 
        0.0711242550022214*Sin[0.017167172970436028*n] + 
        0.0005863132618294766*Sin[0.034334345940872056*n] + 
        0.0013462049383894524*Sin[0.05150151891130809*n]] - 
     1.*Tan[lat]*Tan[0.005782961777094692 - 0.4001419318234436*
         Cos[0.017167172970436028*n] - 0.0060922154967620835*
         Cos[0.034334345940872056*n] - 0.002387468786938206*
         Cos[0.05150151891130809*n] + 0.0711242550022214*
         Sin[0.017167172970436028*n] + 0.0005863132618294766*
         Sin[0.034334345940872056*n] + 0.0013462049383894524*
         Sin[0.05150151891130809*n]]] - 
  2.*ArcCos[-0.00523596383141958*Sec[lat]*Sec[0.005782961777094692 - 
        0.4001419318234436*Cos[0.017167172970436028*n] - 
        0.0060922154967620835*Cos[0.034334345940872056*n] - 
        0.002387468786938206*Cos[0.05150151891130809*n] + 
        0.0711242550022214*Sin[0.017167172970436028*n] + 
        0.0005863132618294766*Sin[0.034334345940872056*n] + 
        0.0013462049383894524*Sin[0.05150151891130809*n]] - 
     1.*Tan[lat]*Tan[0.005782961777094692 - 0.4001419318234436*
         Cos[0.017167172970436028*n] - 0.0060922154967620835*
         Cos[0.034334345940872056*n] - 0.002387468786938206*
         Cos[0.05150151891130809*n] + 0.0711242550022214*
         Sin[0.017167172970436028*n] + 0.0005863132618294766*
         Sin[0.034334345940872056*n] + 0.0013462049383894524*
         Sin[0.05150151891130809*n]]])

(* it is 184 terms long in TreeForm *)

In[158]:= LeafCount[TreeForm[temp1942[n,lat]]]

Out[158]= 184

(*

I'm looking for an approximation that is significantly shorter, and
still reasonably accurate for -Pi/3 <= lat <= Pi/3; the function blows
up as Abs[lat] approaches Pi/2, so the approximation need not be
accurate when Abs[lat] > Pi/3

The approximation should not be a Piecewise function (like FindFormula
tends to yield, and which Interpolate generates by design), and should
be a reasonably easy to calculate "rule of thumb".

The (extremely ugly) work I've done so far is at:

https://github.com/barrycarter/bcapps/blob/master/ASTRO/sun-rise-set-lengths.m

https://github.com/barrycarter/bcapps/blob/master/ASTRO/bc-astro-formulas.m

I'm also looking for general tips on how to approximate a 2
dimensional function that I can evaluate at any point. I've looked
at most of the functions listed on
http://reference.wolfram.com/language/FunctionApproximations/tutorial/FunctionApproximations.html
by fixing one value in my function (usually 'n'), and then trying to
find a pattern to the approximating functions. This does not work
well.

My goal with this specific function is to answer
https://astronomy.stackexchange.com/questions/24304/expression-for-length-of-sunrise-sunset-as-function-of-latitude-and-day-of-year
more succinctly than the current answer, but more generally to provide
good approximations to other (mostly astronomical) formulas that have
no closed form.

*)
$\endgroup$
2
  • $\begingroup$ I'm missing some kind of symmetrie in your function temp1942[n, lat] == temp1942[n+365/366?, lat] . You gave a range 0.5 <= n <= 366.5 but I didn't succeed finding the period. $\endgroup$ Commented Feb 15, 2018 at 9:16
  • $\begingroup$ @UlrichNeumann Note that temp1942[366.5,lat] == temp1942[0.5,lat]; the period is of length 366. The Sin/Cos functions should all have periods that are divisors of 366. $\endgroup$
    – user1722
    Commented Feb 22, 2018 at 3:46

1 Answer 1

0
$\begingroup$

One way to find an approximation is the use of basic functions h which follow from a product approach:

f[\[Alpha]_] := {1, \[Alpha]^2, \[Alpha]^4}
g[n_] := {1, Cos[4 Pi n/366] }
h = Flatten[Outer[Times, f[\[Alpha]], g[n]]] (* basicfunction *)
(*{1, Cos[(2 n \[Pi])/183], \[Alpha]^2, \[Alpha]^2 Cos[(2 n \[Pi])/183], \[Alpha]^4, \[Alpha]^4 Cos[(2 n \[Pi])/183]} *)

The underlying idea is to approximate your function temp1942[n, \[Alpha]]~c.h[n,\[Alpha]] by a linear combination of the basic functios h.

The minimization J=Integrate[(temp1942[n, \[Alpha]]~c.h[n,\[Alpha]])^2,, {n, 0, 366}, {\[Alpha], -Pi/3,Pi/3}]=...=c.M.c-2r.c+const[c] can be solved explicitely for c:

M = NIntegrate[Outer[Times, h, h], {n, 0, 366}, {\[Alpha], -Pi/3,Pi/3},PrecisionGoal -> 5, AccuracyGoal -> 4]
r = NIntegrate[h temp1942[n, \[Alpha]], {n, 0, 366}, {\[Alpha], -Pi/3,Pi/3},PrecisionGoal -> 5, AccuracyGoal -> 4]
c=Inverse[M].r
Plot3D[c.h, {n, 0., 365}, {\[Alpha], -Pi/3, Pi/3},PlotRange -> {0, 600}, MaxRecursion -> 4]

enter image description here

The approximation might be improved using well adapted, additional basic functions.

$\endgroup$
3
  • $\begingroup$ I understand the concept, but didn't 100% follow what you did here, but it works well, thanks! $\endgroup$
    – user1722
    Commented Feb 15, 2018 at 16:31
  • $\begingroup$ The main issue is the minimization of Integrate[(temp1942[n, \[Alpha]]~c.h[n,\[Alpha]])^2,, {n, 0, 366}, {\[Alpha], -Pi/3,Pi/3}](see my edit). This quadratic least square problem can be solved explicitely! Alternatively you could try to solve using NMinimize directly, but the performance would be poor because NMinimize calls NIntegrate some 100times. $\endgroup$ Commented Feb 16, 2018 at 7:10
  • $\begingroup$ Thanks! I eventually figured this out. You're saying: if you approximate f as a1*f1 + a2*f2 + a3*f3 + ..., then the integral of the difference squared, (f - (a1*f1 + a2*f2 + a3*f3 + ..))^2, is just a linear combination of the integrals of f^2, ai * fi * f, and ai * aj * fi * fj. Since the ai are constants, it suffices to find the definite integrals of f^2, fi*f, and fi*fj, which will just be real numbers. Then, you have a linear combination of {1, ai, ai*aj} (including i=j), which can be minimized using least squares. $\endgroup$
    – user1722
    Commented Feb 22, 2018 at 3:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.