       Re: Factor, Expand. Daytime Hours.

• To: mathgroup at smc.vnet.net
• Subject: [mg32320] Re: Factor, Expand. Daytime Hours.
• From: Tom Burton <tburton at cts.com>
• Date: Sat, 12 Jan 2002 21:32:39 -0500 (EST)
• References: <a1p2gh\$6vq\$1@smc.vnet.net>
• Reply-to: tburton at cts.com
• Sender: owner-wri-mathgroup at wolfram.com

```Hi,

On Sat, 12 Jan 2002 10:17:53 +0000 (UTC), in comp.soft-sys.math.mathematica you wrote:

>2.
>Where can I find the function that gives the daytime hours between sunrise
>and sunset as a function of date and latitude?

I might be able to help with your second question. Paste the following expression into a notebook to yield a function sunPath for sun elevation versus latitude, season of year, and time of day. You can easily apply FindRoot to the result to find times of sunrize and sunset.

Most of the notebook is devoted to derivations, examples, illustrative plots, and supporting functions. But you initialize the notebook and follow the usage statement to get going.

Tom Burton

Notebook[{

Cell[CellGroupData[{
Cell[TextData[{
"Prepared by Tom Burton, ",
ButtonBox["Brahea",
ButtonData:>{
URL[ "http://www.brahea.com";], None},
", Inc.\nPrepared for myself\nJuly 2000"
}], "Subsubtitle",
ShowGroupOpenCloseIcon->False],

Cell[CellGroupData[{

Cell["External dependencies", "Section"],

Cell[CellGroupData[{

Cell[TextData[{
"Standard ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" packages"
}], "Subsection"],

Cell[BoxData[
\(Needs["\<Graphics`Colors`\>"]\)], "Input",
InitializationCell->True],

Cell[BoxData[
\(Needs["\<Graphics`Arrow`\>"]\)], "Input",
InitializationCell->True]
}, Open  ]]
}, Closed]],

Cell[CellGroupData[{

Cell["Graphics options", "Section"],

Cell[BoxData[
\(\$DefaultFont = {"\<Arial-Bold\>", 10}\)], "Input",
InitializationCell->True],

Cell[BoxData[
\(Needs["\<Graphics`Graphics`\>"]\)], "Input",
InitializationCell->True],

Cell[BoxData[
\(Needs["\<Graphics`MultipleListPlot`\>"]\)], "Input",
InitializationCell->True],

Cell[BoxData[
\(With[{myLineStyles = {\n (*\
Crv\ \(1 : \ solid\)\ \ \ \ \ \ *) \ {Thickness[
0.006]}, \n (*\ \ \ \ \ 2 : \
long\ dash\ \ *) \ {Thickness[0.008], \
Dashing[{0.03, \ 0.02}]}, \n (*\ \ \ \ \ 3 : \
dot\ dash\ \ \ *) \ {Thickness[
0.008], \ \n\ \ \ \ \ \ \ \ \ \ Dashing[{0.005, \
0.012, \ 0.015, \ 0.012}]}, \n (*\ \ \ \ \ 4 : \
short\ dash\ *) \ {Thickness[0.011], \
Dashing[{0.005, \ 0.015}]}\n (*\ \ \ \ \ 5 : \
solid, \ \(\(etc\)\(.\)\), \
repeats\ ad\ infinitum\ *) \ }}, \n\t
SetOptions[Plot, \n\t\tGridLines \[Rule] Automatic,
PlotRange \[Rule] All, \n\ \ \ \ PlotStyle \[Rule]
myLineStyles]; \n\t
SetOptions[LogLogPlot, \n\t\tGridLines \[Rule] Automatic,
PlotRange \[Rule] All, \n\ \ \ \ PlotStyle \[Rule]
myLineStyles]; \n\t
SetOptions[LogLinearPlot, \n\t\tGridLines \[Rule] Automatic,
PlotRange \[Rule] All, \n\ \ \ \ PlotStyle \[Rule]
myLineStyles]; \n\t
SetOptions[LinearLogPlot, \n\t\tGridLines \[Rule] Automatic,
PlotRange \[Rule] All, \n\ \ \ \ PlotStyle \[Rule]
myLineStyles]; \n\t
SetOptions[MultipleListPlot, \n\t\tGridLines \[Rule] None,
PlotRange \[Rule] All,
PlotStyle \[Rule] myLineStyles];\n]\)], "Input",
InitializationCell->True],

Cell[BoxData[
\(With\ [{ptsz\  = \ 0.015\  (*\ PC : \ 0.015\ for\ screen, \
0.005\ for\ print/export\ *) , \n\ \ \ \ \ \ \ thk\ \  = \
0.003}, \n
SetOptions\ [
ListPlot, \ \ \ \ \ \ \ \ \ \ GridLines\  \[Rule] \
Automatic, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ PlotRange\  \[Rule] \
All, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ PlotStyle\  \[Rule] \ {PointSize[
ptsz], \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[thk]}\ ]\ ; \n
SetOptions\ [LogLinearListPlot, \
GridLines\  \[Rule] \
Automatic, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ PlotRange\  \[Rule] \
All, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ PlotStyle\  \[Rule] \ {PointSize[
ptsz], \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[thk]}\ ]\ ; \n
SetOptions\ [LinearLogListPlot, \
GridLines\  \[Rule] \
Automatic, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ PlotRange\  \[Rule] \
All, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ PlotStyle\  \[Rule] \ {PointSize[
ptsz], \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[thk]}\ ]\ ; \n
SetOptions\ [
LogLogListPlot, \ \ \ \ GridLines\  \[Rule] \
Automatic, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ PlotRange\  \[Rule] \

All, \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ PlotStyle\  \[Rule] \ {PointSize[
ptsz], \n\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Thickness[
thk]}\ ]\ ;\n]\)], "Input",
InitializationCell->True]
}, Closed]],

Cell[CellGroupData[{

Cell["Basis functions", "Section",
PageWidth->WindowWidth],

Cell[TextData[{
"The coordinate system is uniquely is by the directions of two of \
its three basis vectors {",
"x",
",",
"y",
",",
"z",
"}. The basis is a matrix whose ",
"rows",
" are the three basis vectors.The function ",
"basis",
" accepts one or two unnormalized vectors representing ",
"x",
", ",
"y",
", or ",
"z",
". If two of these are provided but are colineaer and therefore \
cannot form a basis, then the global basis ",
Cell[BoxData[
RowBox[{"(", GridBox[{
{"1", "0", "0"},
{"0", "1", "0"},
{"0", "0", "1"}
}], ")"}]]],
" is returned.  There are also special forms defined below."
}], "Text"],

Cell[CellGroupData[{

Cell["Usage and messages", "Subsection",
PageWidth->WindowWidth],

Cell[BoxData[
\(\(basis::"\<usage\>" = "\<Compute a basis (a matrix of three \
orthormal row-vectors)\n
basis_matrix = basis[...] Several input forms:\n
basis[local_x, local_y], makes z_local\n
basis[local_z], axisym beams and iso plates\n
basis[local_x,0,local_z], beams and ortho plates\n
basis[0,local_y,local_z], alternative for beams\n
If needed, orthogonolize local_y w.r.t local_x, local_x w.r.t. \
local_z\>";\)\)], "Input",
PageWidth->WindowWidth,
InitializationCell->True],

Cell[TextData[{
"These routines write the three basis vectors of a coordinates \
system into the ",
StyleBox["rows",
FontSlant->"Italic"],
" of a basis matrix, using a convention compatible with the \
transformations defined in this notebook. Row-matrices are more \
convenient in ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" than column matrices, because ",
StyleBox["Mathematica",
FontSlant->"Italic"],
" stores  matrix by rows. Convential mathematical notation, \
however, is to write these basis vectors into the ",
StyleBox["columns",
FontSlant->"Italic"],
" of the matrix. Beware."
}], "Text",
PageWidth->WindowWidth],

Cell[BoxData[
\(\(basis::parallel = "\<`1`-vector `2` is parallel to `3`-vector \
`4`. The basis cannot be formed.\>";\)\)], "Input",
InitializationCell->True],

Cell[BoxData[
\(\(basis::short = "\<`1`-vector `2` or `3`-vector `4` has zero \
length. The basis cannot be formed.\>";\)\)], "Input",
InitializationCell->True]
}, Open  ]],

Cell[CellGroupData[{

Cell["General basis function", "Subsection",
PageWidth->WindowWidth,
Evaluatable->False],

Cell[TextData[{
"This function returns a ",
StyleBox["row-matrix",
FontSlant->"Italic"],
" of orthonormal vectors based on input vectors local_x and \
local_y. Local_z is constructed."
}], "Text",
PageWidth->WindowWidth,
Evaluatable->False],

Cell[CellGroupData[{

Cell["Code", "Subsubsection",
PageWidth->WindowWidth],

Cell["\<\
v. 030 (7/24/96) Force numerical return (// N)
v. 057 (1/9/98) Abort on exception\
\>", "Text",
PageWidth->WindowWidth,
Evaluatable->False,
FontFamily->"Times New Roman",
FontSize->11,
FontWeight->"Plain",
FontSlant->"Plain",
FontTracking->"Plain",
FontColor->GrayLevel,
Background->GrayLevel,
FontVariations->{"Underline"->False,
"Outline"->False,

Cell[BoxData[
\(basis[a_List, b_List] :=
Module[{anorm, bo, bnorm}, \n\t\tIf[
a . a\ b . b == 0, \n\t\t\tMessage[basis::short, "\<X\>",
a, "\<Y\>", b]; \n\t\t\tAbort\n\t\t]; \n\t\tIf[
a\[Cross]b . a\[Cross]b == 0, \n\t\t\tMessage[
basis::parallel, "\<X\>", a, "\<Y\>", b]; \n\t\t\tAbort[
1]\n\t\t]; \n\t\tanorm = a\/\@\(a . a\);
bo = b - anorm\ anorm . b; bnorm = bo\/\@\(bo . bo\);
N[{anorm, bnorm, anorm\[Cross]bnorm}]]\)], "Input",
PageWidth->WindowWidth,
InitializationCell->True,
FontFamily->"Courier New",
FontWeight->"Bold",
FontColor->GrayLevel,
Background->GrayLevel]
}, Open  ]],

Cell[CellGroupData[{

Cell["Local tests", "Subsubsection",
PageWidth->WindowWidth],

Cell["\<\
basis[{0,0,0},{0,0,1}]
basis[{0,0,1},{0,0,2}]
basis[{0,1,2},{1,2,3}]\
\>", "Input",
PageWidth->WindowWidth,
FontFamily->"Courier New",
FontWeight->"Bold",
FontColor->GrayLevel,
Background->GrayLevel]
}, Closed]]
}, Open  ]],

Cell[CellGroupData[{

Cell["Axisymmetric beam or isotropic plate", "Subsection",
PageWidth->WindowWidth],

Cell["\<\
First case, basis function to orient an isotropic beam. The vector \
connects end 1 to end 2 and is assumed to define the beam's third \
principle direction. The other two are arbitrary. Also used to orient \
a curved connector.\
\>", "Text",
PageWidth->WindowWidth,
Evaluatable->False,
FontFamily->"Times New Roman",
FontSize->11,
FontWeight->"Plain",
FontSlant->"Plain",
FontTracking->"Plain",
FontColor->GrayLevel,
Background->GrayLevel,
FontVariations->{"Underline"->False,
"Outline"->False,

Cell[CellGroupData[{

Cell["Code", "Subsubsection"],

Cell[TextData[{
"v. 030 (7/24/96) Force numerical return (// N)\nv. 057 (1/8/98) \
Abort on an exception\nv. 083 (12/16/98) Change criterion for ",
StyleBox["c", "Input"],
" being close to parallel to ",
StyleBox["b", "Input"]
}], "Text"],

Cell["\<\
basis[c_List] := Module[{cnorm, b={0,1,0}, a},
If[c.c==0,
Message[basis::\"short\",\"Z\",c,\"(NA)\",\"\"];
Abort
];
cnorm=c/Sqrt[c.c];
(* Assign the second basis vector to {0,1,0} unless that is too \
close
to being parallel to c *)
If[ Abs[cnorm.b] > 0.7, b = {1,0,0} ];
b = b - cnorm cnorm.b;        (* orthogonalize b w.r.t. c *)
b = b/Sqrt[b.b];
a =  b\[Cross]cnorm;
{a, b, cnorm} // N
]\
\>", "Input",
PageWidth->WindowWidth,
InitializationCell->True,
FontFamily->"Courier New",
FontWeight->"Bold",
FontColor->GrayLevel,
Background->GrayLevel]
}, Open  ]],

Cell[CellGroupData[{

Cell["Local tests", "Subsubsection",
PageWidth->WindowWidth],

Cell["\<\
basis[{0,0,0}]
basis[{0,0,10}]
basis[{0,10,0}]
basis[{0,-1.5,0}]
basis[{3,2,1}]\
\>", "Input",
PageWidth->WindowWidth,
FontFamily->"Courier New",
FontSize->10,
FontWeight->"Bold",
FontColor->GrayLevel,
Background->GrayLevel],

Cell[BoxData[
\(monitorReset[]\)], "Input"]
}, Open  ]]
}, Closed]],

Cell[CellGroupData[{

Cell["Orthotropic beam basis functions", "Subsection",
PageWidth->WindowWidth],

Cell["\<\
Second case, orthotropic beam. This function is also used by the \
orthotropic plate, specified by first material direction and plate \
normal direction.\
\>", "Text",
PageWidth->WindowWidth],

Cell[TextData[{
"The third principal direction of the beam is assumed to align with \
vector ",
StyleBox["c",
FontSlant->"Italic"],
". Vector a points in the first principal direction. It will be \
orthogonalized against ",
StyleBox["c",
FontSlant->"Italic"],
" just be sure."
}], "Text",
PageWidth->WindowWidth],

Cell[CellGroupData[{

Cell["Code", "Subsubsection",
PageWidth->WindowWidth],

Cell["\<\
v. 030 (7/24/96) Force numerical return (// N)
v. 055 (9/21/97) Orthogonalize a w.r.t. c instead of c w.r.t a
v. 057 (1/8/98) Abort on exceptions\
\>", "Text",
PageWidth->WindowWidth],

Cell[BoxData[
\(basis[a_List, 0, c_List] :=
Module[{anorm, b, ao, cnorm}, \n\t\tIf[
a . a\ c . c == 0, \n\t\t\t\ \ Message[
basis::"\<short\>", "\<X\>", a, "\<Z\>",
c]; \n\t\t\tAbort\n\t\t]; \n\t\tIf[
a\[Cross]c == {0, 0, 0}, \n\t\t\t\ \ Message[
basis::"\<parallel\>", "\<X\>", a, "\<Z\>",
c]; \n\t\t\tAbort\n\t\t]; \n\t\tcnorm =
c\/\@\(c . c\); ao = a - cnorm\ cnorm . a;
anorm = ao\/\@\(ao . ao\); b = cnorm\[Cross]anorm;
N[{anorm, b, cnorm}]]\)], "Input",
PageWidth->WindowWidth,
InitializationCell->True],

Cell["\<\
v. 057 (12/24/97) new function
v. 057 (1/8/98) Abort on exceptions\
\>", "Text"],

Cell[BoxData[
\(basis[0, b_List, c_List] :=
Module[{a, bnorm, bo, cnorm}, \n\t\tIf[
b . b\ c . c == 0, \n\t\t\tMessage[
basis::"\<short\>", "\<Y\>", b, "\<Z\>",
c]; \n\t\t\tAbort\n\t\t]; \n\t\tIf[
b\[Cross]c == {0, 0, 0}, \n\t\t\tMessage[
basis::"\<parallel\>", "\<Y\>", b, "\<Z\>",
c]; \n\t\t\tAbort\n\t\t]; \n\t\tcnorm =
c\/\@\(c . c\); bo = b - cnorm\ cnorm . b;
bnorm = bo\/\@\(bo . bo\); a = bnorm\[Cross]cnorm;
N[{a, bnorm, cnorm}]]\)], "Input",
PageWidth->WindowWidth,
InitializationCell->True]
}, Open  ]],

Cell[CellGroupData[{

Cell["Local tests", "Subsubsection",
PageWidth->WindowWidth],

Cell["\<\
basis[{0,0,3},0,{0,0,0}]
basis[{0,0,0},0,{0,0,4}]
basis[{0,0,3},0,{0,4,100}]
basis[{3,0,0},0,{3,0,1}]
basis[{0,0,3},0,{100,4,0}]\
\>", "Input",
PageWidth->WindowWidth],

Cell["\<\
basis[0,{0,0,3},{0,0,2}]
basis[0,{0,0,0},{0,0,4}]
basis[0,{0,0,3},{0,4,100}]
basis[0,{3,0,0},{3,0,1}]
basis[0,{0,0,3},{100,4,0}]\
\>", "Input",
PageWidth->WindowWidth]
}, Closed]]
}, Open  ]],

Cell[CellGroupData[{

Cell["Basis of a line segment", "Subsection"],

Cell[TextData[{
"The general shape of a line junction is a piecewise-straight line \
in 3-space. The local ",
StyleBox["x",
FontSlant->"Italic"],
"-axes of one straight segment of a line junction and all of its \
connections are aligned with each other and with the segment itself. \
The ",
StyleBox["AutoSEA",
FontSlant->"Italic"],
" engine characterizes the orientation of a segment of a line \
connection in terms of the angle \[Beta] that the local ",
StyleBox["y",
FontSlant->"Italic"],
"-axis of the connection's segment makes with that of the \
junction's segment. Then the orientation of the connection as a whole \
is the average of \[Beta] over the connection, weighted by arclength. \
(Look ",
ButtonBox["here",
ButtonData:>{"Globals.nb", "cxn basis"},
" for futher details.) This approach requires that the local basis \
of every segment of every line junction be uniquely defined. This \
section lays out one such definition. Other possible definitions \
would yield the same results, provided that they are also unique."
}], "Text"],

Cell[TextData[{
"The local row-basis {",
StyleBox["X",
FontWeight->"Bold"],
", ",
StyleBox["Y",
FontWeight->"Bold"],
", ",
StyleBox["Z",
FontWeight->"Bold"],
"} of a straight segment of line junction will be unique if it is \
defined unambiquously in relation to a single global row-basis \
{Xg,Yg, Zg}. Without loss of generality, we can choose {Xg,Yg, Zg}=",
Cell[BoxData[
RowBox[{"(", GridBox[{
{"1", "0", "0"},
{"0", "1", "0"},
{"0", "0", "1"}
}], ")"}]]],
". Note that ",
StyleBox["X",
FontWeight->"Bold"],
" is fixed; only ",
StyleBox["Y",
FontWeight->"Bold"],
" and ",
StyleBox["Z",
FontWeight->"Bold"],
" need be determined. Then here is the unique definition that we \
have chosen:"
}], "Text"],

Cell[CellGroupData[{

Cell["Usage", "Subsubsection"],

Cell[BoxData[
\(\(basisLine::"\<usage\>" = "\<basisLine[X_vector] \
\[LongRightArrow] {X,Y,Z}, a unique basis for a line segment parallel \
to X.\>";\)\)], "Input",
InitializationCell->True]
}, Open  ]],

Cell[CellGroupData[{

Cell["Code", "Subsubsection"],

Cell["v. 057 (12/13/97) new function", "Text"],

Cell[BoxData[
RowBox[{
RowBox[{"basisLine", "=",
RowBox[{"Compile", "[",
RowBox[{\({{a, _Real, 1}}\), ",",
RowBox[{"Module", "[",

RowBox[{\({X = a/\@\(a . a\), Xg, Yg, \ Zg, b, c, bb,
cc}\), ",", "\n", "\t\t",
RowBox[{
RowBox[{\({Xg, Yg, \ Zg}\), "=",
RowBox[{"(", GridBox[{
{"1", "0", "0"},
{"0", "1", "0"},
{"0", "0", "1"}
}], ")"}]}], ";", \(b = Zg\[Cross]X\),
";", \(c = X\[Cross]Yg\), ";", \(bb = b . b\),
";", \(cc = c . c\), ";", "\n",
"\t\t", \(If[bb > cc, \n\t\t\tY = b/\@bb;
Z = X\[Cross]Y, \n\t\t\tZ = c/\@cc;
Y = Z\[Cross]X\n\t\t]\), ";", "\n",
"\t\t", \({X, Y, Z}\)}]}], "\n", "\t", "]"}]}],
"]"}]}], ";"}]], "Input",
InitializationCell->True],

Cell["\<\
Here are four pictures showing the bases generated for sequences of \
line arranged in circles of four different orientations. Bases are \
represented by trials with {black, red, blue} representing {X,Y,Z}. \
Note that, in each case, the black lines (X vectors) form a circle.\
\>", "Text"],

Cell[BoxData[
RowBox[{
RowBox[{"Show", "[",
RowBox[{
RowBox[{"GraphicsArray", "[", "\n", "\t\t",
RowBox[{"Apply", "[", "\n", "\t\t\t",

RowBox[{\(With[{pts =
Partition[
Table[#1\ Cos[\[Theta]\ Degree] + #2\ Sin[\
\[Theta]\ Degree], {\[Theta], 0. , 360. , 20. }], 2,

1]}, \n\t\t\t\tWith[{d = \(0.6 \@\( # . #\) \
&\)[pts\[LeftDoubleBracket]1, 2\[RightDoubleBracket] -
pts\[LeftDoubleBracket]1,
1\[RightDoubleBracket]]}, \n\t\t\t\t\t\
Graphics3D[Function[{pq}, {p, q} =
pq; \n\t\t\t\t\t\t\tTranspose[{{black,
red,
blue}, \(Line[{p, p + d\ #}] &\) /@
basisLine[q - p]}]] /@ pts]]] &\),
",", "\n", "\t\t\t",
RowBox[{"(", GridBox[{
{\({{1, 0, 0}, {0, 1, 0}}\), \({{0, 1, 0}, {0,
0, 1}}\)},
{\({{0, 0, 1}, {1, 0, 0}}\), \({{1, 1, 1}, {1,
0, \(-1\)}}\)}
}], ")"}], ",", "\n", "\t\t\t", \({2}\)}],
"]"}], "]"}], ",", "\n",
"\t", \(ImageSize \[Rule] {500, 400}\)}], "]"}],
";"}]], "Input"],

Cell["\<\
In the left-hand pair of diagrams, the triad moves smoothly around \
the circle. But in the right-hand pair, the trial rotates abruptly. \
This kind of abrupt rotation is hard to or impossible to avoid in \
general, but is acceptable if it happens consistently and repeatably.\
\
\>", "Text"]
}, Closed]]
}, Closed]]
}, Closed]],

Cell[CellGroupData[{

Cell["Sun path", "Section"],

Cell["\<\
The sunPath function returns the azimuth \[Theta] and altitude \[Phi] \
of the sun versus latitude, and season, and sun time.\
\>", "Text"],

Cell[BoxData[
\(\(sunPath::usage = "\<sunPath[latitude, season, suntime] \
\[LongRightArrow] {azimuth, altitude} of the sun. All parameters are \
angles in radians. Season is zero and \[Pi] at the Northern winter \
and summer solstices. Suntime is zero at midnight at zero longitude \
and \[Pi] at high noon there. The earth inclination of 23.45\[Degree] \
is assumed.\>";\)\)], "Input",
InitializationCell->True],

Cell[CellGroupData[{

Cell["Development", "Subsection"],

Cell[TextData[{
"Let the orbit of the earth be in the global ",
StyleBox["XY",
FontSlant->"Italic"],
" plane. Orient this plane such that the earth's axis of rotation \
is in the ",
StyleBox["XZ",
FontSlant->"Italic"],
" plane. Call by the name \[Psi] this angle of inclination of the \
earth's axis relative to ",
StyleBox["Z.",
FontSlant->"Italic"],
"  Let the angle of the earth around the sun\[LongDash]its position \
in orbit\[LongDash]with respect to the global ",
StyleBox["X",
FontSlant->"Italic"],
" axis be \[Alpha]. Call the vector from sun to earth and its \
magnitude ",
Cell[BoxData[
" and ",
Cell[BoxData[
". Then ",
Cell[BoxData[
R[\[Alpha]] {Cos[\[Alpha]], Sin[\[Alpha]], 0}\)]],
"."
}], "Text"],

Cell[BoxData[
\(\(Show[
Graphics[{Yellow, Disk[{0, 0}, 10], Black,
Line[\({100, 0} + # { .3, 1} &\) /@ {\(-6\), 6}],
Text["\<z\>", {100, 0} + 9 { .3, 1}],
Line[\({100, 0} + # {1, \(-0.3\)} &\) /@ {\(-1\), 6}],
Text["\<x\>", {100, 0} + 8 {1, \(- .3\)}], Blue,
Disk[{100, 0}, 2], Black, Thickness[0.005],
Arrow[{0, \(-15\)}, {100, \(-15\)}],
Text[R\&\[RightVector], {50, \(-15\)}, {0, \(-1\)}]\
\[IndentingNewLine]}, AspectRatio \[Rule] Automatic,
PlotRange \[Rule] {{\(-10\), 110}, {\(-20\), 15}},
Axes \[Rule] True, Ticks \[Rule] None,
AxesLabel \[Rule] {"\<X\>", "\<Z\>"},
PlotLabel \[Rule] "\<Midnight near the Winter solstice\>",
ImageSize \[Rule] 450]];\)\)], "Input",
CellTags->"hide"],

Cell[TextData[{
"The global Cartesian coordinate system ",
StyleBox["XYZ",
FontSlant->"Italic"],
" of the solar system and the system of the earth ",
StyleBox["xyz",
FontSlant->"Italic"],
" at m",
"idnight at 0\[Degree] longitude near the Winter solstice of the \
Northern Hemisphere."
}], "Text",
CellMargins->{{72, 80}, {Inherited, Inherited}},
FontFamily->"Arial",
FontSize->10,
FontWeight->"Bold"],

Cell[TextData[{
"Let a coordinate system embedded in the earth be called earth with \
",
StyleBox["x",
FontSlant->"Italic"],
" in the ",
StyleBox["RZ",
FontSlant->"Italic"],
" plane away from the sun and ",
StyleBox["z",
FontSlant->"Italic"],
" along the axis \[PlusMinus]90\[Degree] latitude. This basis may \
be expressed in global coordinates as follows"
}], "Text"],

Cell[BoxData[
\(\(\(\[IndentingNewLine]\)\(Map[
Composition[Rationalize, FullSimplify],
basis[{Cos[\[Alpha]], Sin[\[Alpha]], 0},
0, {Sin[\[Psi]], 0,
Cos[\[Psi]]}], \[IndentingNewLine]\(-1\)]\)\)\)], "Input"],

Cell[TextData[{
"The inclination angle is well known to be ",
Cell[BoxData[
\(TraditionalForm\`\[Psi] = 23.45  \[Degree]\)]],
". so we have"
}], "Text"],

Cell[BoxData[
\(earth =
With[{\[Psi] = 23.45  \[Degree]}, \[IndentingNewLine]Map[
Composition[Rationalize, FullSimplify],
basis[{Cos[\[Alpha]], Sin[\[Alpha]], 0},
0, {Sin[\[Psi]], 0,
Cos[\[Psi]]}], \[IndentingNewLine]\(-1\)]\
\[IndentingNewLine]]\)], "Input"],

Cell[TextData[{
"Let the rotation of the earth be described in terms of an angle \
\[Beta] between longitude 0\[Degree] and the global ",
StyleBox["RZ",
FontSlant->"Italic"],
" plane. Now erect a local coordinate system at 20\[Degree] North \
latitude and 155\[Degree] West longitude (close to Hakalau) with ",
StyleBox["x",
FontSlant->"Italic"],
" pointing East, ",
StyleBox["y",
FontSlant->"Italic"],
" pointing North, and ",
StyleBox["z",
FontSlant->"Italic"],
" pointing up. Its basis in coordinates of the earth is"
}], "Text"],

Cell[BoxData[
\(\(local =
With[{\[CapitalTheta] = \(-155.0\) \[Degree] + \[Beta], \
\[CapitalPhi] = \(+20.0\) \[Degree]}, \
\[IndentingNewLine]basis[{\(-Sin[\[CapitalTheta]]\),
Cos[\[CapitalTheta]], 0},
0, {Cos[\[CapitalTheta]] Cos[\[CapitalPhi]],
Sin[\[CapitalTheta]] Cos[\[CapitalPhi]],
Sin[\[CapitalPhi]]}]\[IndentingNewLine]];\)\)], \
"Input"],

Cell[BoxData[
\(FullSimplify@
Map[Rationalize,
With[{\[CapitalTheta] = \[Beta]}, \
\[IndentingNewLine]basis[{\(-Sin[\[CapitalTheta]]\),
Cos[\[CapitalTheta]], 0},
0, {Cos[\[CapitalTheta]] Cos[\[CapitalPhi]],
Sin[\[CapitalTheta]] Cos[\[CapitalPhi]],
Sin[\[CapitalPhi]]}]\[IndentingNewLine]], \(-1\)]\)], \
"Input"],

Cell[TextData[{
"It is midnight, sun time, when ",
Cell[BoxData[
"."
}], "Text"],

Cell["\<\
For reference, the prematrix rotating a vector from the global \
coordinate system to another system is\
\>", "Text"],

Cell[BoxData[
\(Clear[a, b, c]\)], "Input"],

Cell[BoxData[
RowBox[{"G3s", "[",
RowBox[{\(IdentityMatrix\), ",",
RowBox[{"(", GridBox[{
{"a", "b", "c"},
{"d", "e", "f"},
{"g", "h", "i"}
}], ")"}]}], "]"}]], "Input"],

Cell["Thus, the prematrix converting from global to local is", \
"Text"],

Cell[BoxData[
\(Rationalize[local . earth]\)], "Input"],

Cell["\<\
Recall that the unit vector from the earth to the sun in global \
coordinates is\
\>", "Text"],

Cell[BoxData[
\({Cos[\[Alpha] + \[Pi]], Sin[\[Alpha] + \[Pi]], 0}\)], "Input"],

Cell["Therefore this vector in local coordinates is", "Text"],

Cell[BoxData[
\({a, b, c} =
Together[
Rationalize[
local . earth] . {\(-Cos[\[Alpha]]\), \(-Sin[\[Alpha]]\),
0}]\)], "Input"],

Cell[BoxData[
\(\[Theta] = ArcTan[a, b]\)], "Input"],

Cell["and the latitude of this vector is", "Text"],

Cell[BoxData[
\(\[Phi] = ArcSin[c]\)], "Input"],

Cell["At high noon on the winter solstice,", "Text"],

Cell[BoxData[
\(Block[{\[Alpha] =
0, \[Beta] = \[Pi] +
155  \[Degree]}, \[IndentingNewLine]{{a, b,
c}, {\[Theta], \[Phi]}}\[IndentingNewLine]]\)], "Input"],

Cell["At high noon on summer solstice,", "Text"],

Cell[BoxData[
\(Block[{\[Alpha] = \[Pi], \[Beta] = \[Pi] +
155  \[Degree]}, \[IndentingNewLine]{{a, b,
c}, {\[Theta], \[Phi]}}\[IndentingNewLine]]\)], "Input"],

Cell[BoxData[
\(\(Block[{\[Alpha] = 0, \[Beta]}, \[IndentingNewLine]Show[
GraphicsArray[{\[IndentingNewLine]Plot[\[Beta] =
b\ \[Degree] + 155  \[Degree]; \[Theta]/Degree, {b,
90, 270},
AxesLabel \[Rule] {"\<\[Beta]\>", "\<\[Theta]\>"},
DisplayFunction \[Rule] Identity],
Plot[\[Beta] = b\ \[Degree] + 155  \[Degree]; \[Phi]/
Degree, {b, 90, 270},
AxesLabel \[Rule] {"\<\[Beta]\>", "\<\[Phi]\>"},
DisplayFunction \[Rule]
Identity]\[IndentingNewLine]},
PlotLabel \[Rule] "\<Winter solstice\>",
ImageSize \[Rule]
500]]\[IndentingNewLine]];\)\)], "Input"],

Cell[BoxData[
\(\(Block[{\[Alpha] = \[Pi], \[Beta]}, \[IndentingNewLine]Show[
GraphicsArray[{\[IndentingNewLine]Plot[\[Beta] =
b\ \[Degree] + 155  \[Degree]; \[Theta]/Degree, {b,
90, 270},
AxesLabel \[Rule] {"\<\[Beta]\>", "\<\[Theta]\>"},
DisplayFunction \[Rule] Identity],
Plot[\[Beta] = b\ \[Degree] + 155  \[Degree]; \[Phi]/
Degree, {b, 90, 270},
AxesLabel \[Rule] {"\<\[Beta]\>", "\<\[Phi]\>"},
DisplayFunction \[Rule]
Identity]\[IndentingNewLine]},
PlotLabel \[Rule] "\<Summer solstice\>",
ImageSize \[Rule]
500]]\[IndentingNewLine]];\)\)], "Input"],

Cell["Clean up.", "Text"],

Cell[BoxData[
\(Clear[a, b, c, local, earth, \[Theta], \[Phi]]\)], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Code", "Subsection"],

Cell["\<\
Set to zero the local longitude, which is cancelling out in the \
development.\
\>", "Text"],

Cell[BoxData[
\(sunPath[\[CapitalPhi]_, \[Alpha]_, \[Beta]_] \
:= \[IndentingNewLine]Module[{earth, local, a, b,

c, \[Theta], \[Phi]}, \[IndentingNewLine]earth = \
{{\(0.8416368868403995`\ Cos[\[Alpha]]\)\/\@\(\(\(0.9208184434201998`\
\)\(\[InvisibleSpace]\)\) - 0.0791815565798002`\ Cos[2\ \[Alpha]]\),
Sin[\[Alpha]]\/\@\(\(\(0.9208184434201998`\)\(\
\[InvisibleSpace]\)\) - 0.0791815565798002`\ Cos[2\ \[Alpha]]\), \
\(-\(\(0.36508113831037636`\ \
Cos[\[Alpha]]\)\/\@\(\(\(0.9208184434201998`\)\(\[InvisibleSpace]\)\) \
- 0.0791815565798002`\ Cos[2\ \[Alpha]]\)\)\)}, \
{\(-\(\(0.9174076993574884`\ \
Sin[\[Alpha]]\)\/\@\(\(\(0.9208184434201998`\)\(\[InvisibleSpace]\)\) \
- 0.0791815565798002`\ Cos[2\ \[Alpha]]\)\)\), \(0.9174076993574884`\ \
Cos[\[Alpha]]\)\/\@\(\(\(0.9208184434201998`\)\(\[InvisibleSpace]\)\) \
- 0.0791815565798002`\ Cos[2\ \[Alpha]]\), \(0.3979486313076105`\ \
Sin[\[Alpha]]\)\/\@\(\(\(0.9208184434201998`\)\(\[InvisibleSpace]\)\) \
- 0.0791815565798002`\ Cos[2\ \[Alpha]]\)}, {0.3979486313076105`, 0,
0.9174076993574884`}}; \[IndentingNewLine]local = \
{{\(-Sin[\[Beta]]\), Cos[\[Beta]],
0}, {\(-Cos[\[Beta]]\)\ Sin[\[CapitalPhi]], \(-Sin[\
\[Beta]]\)\ Sin[\[CapitalPhi]],
Cos[\[CapitalPhi]]}, {Cos[\[Beta]]\ Cos[\[CapitalPhi]],
Cos[\[CapitalPhi]]\ Sin[\[Beta]],
Sin[\[CapitalPhi]]}}; \[IndentingNewLine]{a, b, c} =
local .
earth . {\(-Cos[\[Alpha]]\), \(-Sin[\[Alpha]]\),
0}; \[IndentingNewLine]\[Theta] =
ArcTan[a, b]; \[Phi] =
ArcSin[c]; \[IndentingNewLine]{\[Theta], \[Phi]}\
\[IndentingNewLine]]\)], "Input",
InitializationCell->True]
}, Closed]],

Cell[CellGroupData[{

Cell["Examples", "Subsection"],

Cell[BoxData[
\(sunPath[20  \[Degree], 1, 2]\)], "Input"],

Cell[CellGroupData[{

Cell["Analytical formula", "Subsubsection"],

Cell[BoxData[
\(sunPath[\[CapitalPhi], \[Alpha], \[Beta]]\)], "Input"]
}, Closed]],

Cell[CellGroupData[{

Cell["Sun-path plots", "Subsubsection"],

Cell["\<\
Alter the latitue \[CapitalPhi] in the first argument of Module below \
to suit. The code will not work at latidudes of perpetual winter \
night.\
\>", "Text"],

Cell[BoxData[
\(\(\(\[IndentingNewLine]\)\(figure["\<sun path Hakalau\>"] =
Module[{\[Theta], \[Phi], \[CapitalPhi] =
19.9  \[Degree]}, \[IndentingNewLine]DisplayTogether[\
\[IndentingNewLine]ParametricPlot[
Evaluate[
Table[{\[Theta], \[Phi]} =
sunPath[\[CapitalPhi], \[Alpha], \[Beta]]; \(\((\
\[Pi]/2 - \[Phi])\)\/Degree\) {Cos[\[Theta]],
Sin[\[Theta]]}, {\[Alpha], 0, \ \[Pi]\ ,
2  \[Pi]/12}]], {\[Beta], \[Pi]/2, 3  \[Pi]/2},
Axes \[Rule] False,
PlotRange \[Rule] {{\(-115\), 115}, {\(-110\), 110}},
AspectRatio \[Rule] Automatic,
Epilog \[Rule] {Blue,
Table[Circle[{0, 0}, r], {r, 10, 80, 10}], Black,
Thickness[0.005], Circle[{0, 0}, 90],
Circle[{0, 0},
93], \(Text[\((90 - #)\) \[Degree], {0, #},
Background \[Rule] White] &\) /@
Range[10, 80,
10], \(Text[
Mod[360 + 90 - #, 360, 1] \[Degree],
103 {Cos[#  \[Degree]],
Sin[#  \[Degree]]}] &\) /@ \((Range[30,
360, 30])\),
Text["\<summer\>", {\(-30\), 12}],
Text["\<solstice\>", {30, 12}],
Text["\<winter\>", {\(-30\), \(-45\)}],
Text["\<solstice\>", {30, \(-45\)}], \
\({Text["\<evening\>", {\(-65\), #}],
Text["\<morning\>", {65, #}]} &\) /@ {35, \
\(-45\)}, Text["\<West\>", {\(-109\), 0}, {1, 0}, {0, 1}],
Text["\<East\>", {109, 0}, {\(-1\),
0}, {0, \(-1\)}]},
PlotLabel \[Rule] "\<Sun Path, \>" <>
ToString[\[CapitalPhi]/
Degree] <> "\<\[Degree] North Latitude\>",
ImageSize \[Rule]
450], \[IndentingNewLine]ParametricPlot[
Evaluate[
Table[{\[Theta], \[Phi]} =
sunPath[\[CapitalPhi], \[Alpha], \[Beta]]; \(\((\
\[Pi]/2 - \[Phi])\)\/Degree\) {Cos[\[Theta]],
Sin[\[Theta]]}, {\[Beta], \[Pi]/2,
3  \[Pi]/2, \[Pi]/12}]], {\[Alpha],
0, \[Pi]}], \[IndentingNewLine]ParametricPlot[
Evaluate[
Table[r {Cos[s], Sin[s]}, {s, 0,
2  \[Pi] - \[Pi]/6, \[Pi]/6}]], {r, 90, 95},
PlotStyle \[Rule]
Thickness[0.005]], \[IndentingNewLine]ParametricPlot[
Evaluate[
Table[
r {Cos[s], Sin[s]}, {s, 0,
2  \[Pi] - \[Pi]/18, \[Pi]/18}]], {r, 90, 93},
PlotStyle \[Rule]
Thickness[
0.005]]\[IndentingNewLine]]\[IndentingNewLine]];\)\)\
\)], "Input"],

Cell[BoxData[
\(\(\(\[IndentingNewLine]\)\(figure["\<sun path Santa Fe\>"] =
Module[{\[Theta], \[Phi], \[CapitalPhi] =
35.7  \[Degree]}, \[IndentingNewLine]DisplayTogether[\
\[IndentingNewLine]ParametricPlot[
Evaluate[
Table[{\[Theta], \[Phi]} =
sunPath[\[CapitalPhi], \[Alpha], \[Beta]]; \(\((\
\[Pi]/2 - \[Phi])\)\/Degree\) {Cos[\[Theta]],
Sin[\[Theta]]}, {\[Alpha], 0, \ \[Pi]\ ,
2  \[Pi]/12}]], {\[Beta], \[Pi]/2, 3  \[Pi]/2},
Axes \[Rule] False,
PlotRange \[Rule] {{\(-115\), 115}, {\(-110\), 110}},
AspectRatio \[Rule] Automatic,
Epilog \[Rule] {Blue,
Table[Circle[{0, 0}, r], {r, 10, 80, 10}],
Thickness[0.005], GrayLevel, Circle[{0, 0}, 90],
Circle[{0, 0},
93], \(Text[\((90 - #)\) \[Degree], {0, #},
Background \[Rule] White] &\) /@
Range[10, 80,
10], \(Text[
Mod[360 + 90 - #, 360, 1] \[Degree],
103 {Cos[#  \[Degree]],
Sin[#  \[Degree]]}] &\) /@ \((Range[30,
360, 30])\),
Text["\<summer solstice\>", {0, \(-05\)}],
Text["\<winter solstice\>", {0, \(-65\)}], \
\({Text["\<evening\>", {\(-55\), #}],
Text["\<morning\>", {55, #}]} &\) /@ {25, \
\(-60\)}, Text["\<West\>", {\(-109\), 0}, {1, 0}, {0, 1}],
Text["\<East\>", {109, 0}, {\(-1\),
0}, {0, \(-1\)}]},
PlotLabel \[Rule] "\<Sun Path, \>" <>
ToString[\[CapitalPhi]/
Degree] <> "\<\[Degree] North Latitude\>",
ImageSize \[Rule]
450], \[IndentingNewLine]ParametricPlot[
Evaluate[
Table[{\[Theta], \[Phi]} =
sunPath[\[CapitalPhi], \[Alpha], \[Beta]]; \(\((\
\[Pi]/2 - \[Phi])\)\/Degree\) {Cos[\[Theta]],
Sin[\[Theta]]}, {\[Beta], \[Pi]/2,
3  \[Pi]/2, \[Pi]/12}]], {\[Alpha],
0, \[Pi]}], \[IndentingNewLine]ParametricPlot[
Evaluate[
Table[r {Cos[s], Sin[s]}, {s, 0,
2  \[Pi] - \[Pi]/6, \[Pi]/6}]], {r, 90, 95},
PlotStyle \[Rule]
Thickness[0.005]], \[IndentingNewLine]ParametricPlot[
Evaluate[
Table[r {Cos[s], Sin[s]}, {s, 0,
2  \[Pi] - \[Pi]/18, \[Pi]/18}]], {r, 90, 93},
PlotStyle \[Rule]
Thickness[
0.005]]\[IndentingNewLine]]\[IndentingNewLine]];\)\)\
\)], "Input"]
}, Closed]]
}, Closed]]
}, Open  ]]
}, Open  ]]
},
FrontEndVersion->"4.1 for Microsoft Windows",
ScreenRectangle->{{0, 1280}, {0, 979}},
AutoGeneratedPackage->Automatic,
WindowSize->{975, 740},
WindowMargins->{{Automatic, 34}, {3, Automatic}}
]

```

• Prev by Date: Re: Graphics with column
• Next by Date: Re: Wrestling with Mathematica on Trig Simplification
• Previous by thread: Re: Factor, Expand. Daytime Hours.
• Next by thread: Re: Factor, Expand. Daytime Hours.