How to approximate specific 2 dimensional function + general tips?

Question

(* 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.

*)

Answer 1

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]

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