Fifty Shades of J/Chapter 48

Table of Contents ... Glossary ... Previous Chapter ... Next Chapter

Heavens above!

Principal Topics

o. (circle functions), \. (outfix), spherical trigonometry, similitude, rotations, enlargements, determinant, minors, cofactors, cross product, identity matrix, direction cosines, altitude, azimuth, declination, right ascension, celestial sphere, hour angle, celestial meridian, transit.

Similitude

A general objective of J-ottings has been to draw attention to the considerable number of mathematical or mathematical type routines which are built into J primitives thereby leading to significant reductions of programming effort. One such feature is the versatility of j. which, although primarily a complex number constructor, is adaptable to other circumstances in which objects are defined by pairs of numbers, for example betting odds (see J is complex? You bet!).

A further example concerns transformations of a 2-dimensional plane by means of matrices of the form ${\textstyle M={\begin{pmatrix}a&-b\\b&a\end{pmatrix}}}$ where a and are b real numbers. A transformation such as ${\textstyle M{\begin{pmatrix}x\\y\end{pmatrix}}}$ is called a similitude, that is a transformation which results in the combination of an anti-clockwise rotation about the origin and an enlargement of the objects described by the coordinates. (For a clockwise rotation exchange b and –b ). Given det=:-/ .* standing for determinant, det M is a² + b² and its square root is the enlargement E. The rotation component is represented by M divided by E, resulting in a matrix of the form ${\textstyle M={\begin{pmatrix}\cos t&-\sin t\\\sin t&\cos t\end{pmatrix}}}$  where t is the anti-clockwise angle of rotation.

Now although a similitude could be applied to a triangle whose points are, say, (0,0), (2,1) and (0,1) by

   ]M=.2 2$2 _3 3 2        NB. similitude matrix
2 _3
3  2
    ]tri=.2 3$0 2 0 0 1 1
0 2 0
0 1 1
   M +/ .*tri
0 1 _3
0 8  2

to give the transformed triangle (0,0), (1,8), (-3,2), clearly M is defined by the number pair (a,b) and so can be represented compactly as a J pair, in which case enlargement E and rotation t are given by:

   10 o. 2j3
3.60555               NB. enlargement=sqrt of 2^2 + 3^2
   (% 10&o.)2j3
0.5547j0.83205        NB. (cos x)j(sin x) where tan x=3%2

Instead of using a matrix inner product to transform points, simple multiplication is all that is required, so that the previous triangle transformation is given by:

   j./tri
0 2j1 0j1
   2j3 * 0 2j1 0j1    NB. triangle transformation
0 1j8 _3j2

Also since multiplication is commutative, a product such as 2j3*2j1 has two geometric interpretations, viz. the similitude 2j3 transforms the point (2,1) to the point (1,8) and the similitude 2j1 transforms the point (2,3) to (1,8).

For rotation without enlargement the transformed coordinates of the point (x,y) in the new frame of reference are

$${\begin{pmatrix}\cos t&-\sin t\\\sin t&\cos t\end{pmatrix}}{\begin{pmatrix}x\\y\end{pmatrix}}={\begin{pmatrix}x\cos t-y\sin t\\x\sin t+y\cos t\end{pmatrix}}$$

which can be confirmed by elementary trigonometry:

The components of displacement from (x, y) are thus

$${\begin{pmatrix}x\cos t-y\sin t-x\\x\sin t+y\cos t-y\end{pmatrix}}$$

  or  

$${\begin{pmatrix}x(\cos t-1)-y\sin t\\x\sin t+(y\cos t-1)\end{pmatrix}}$$

which can be written as $-(1-\cos t){\begin{pmatrix}x\\y\end{pmatrix}}-\sin t{\begin{pmatrix}y\\-x\end{pmatrix}}$. The object of this rearrangement will become apparent shortly.

Rotations in three dimensions

Here the geometry is more complicated and j. no longer helps. Take those rotations in which any line through the origin may be chosen as axis. Define such a rotation by any point on it other than the origin, and normalise this so that the defining point lies on the unit sphere. The results of this normalisation are the direction cosines of the axis of rotation, that is the cosines of the angles which this axis makes separately with each of the coordinate axes:

   dircos=:% %:@(+/@:*:)        NB. direction cosines, i.e. normalise rotation axis to unit vector
   dircos 3 4 5
0.424264 0.565685 0.707107

A necessary preliminary is to obtain the cross-product of the rotation vector and a point to be rotated. To my knowledge there is no primitive which delivers cross-products directly, however it is a reasonably straightforward to write a verb xp. First stitch 3 4 5 (defining the axis) to 1 2 3, a point to be rotated, and use the ‘all but one’ technique described in J-ottings 52 to obtain the submatrices obtained by progressively eliminating one row at a time (1 ]\.). The determinant of each of these 2x2 matrices is required with a suitable adjustment for alternating signs leading to:

   xp=:4 : '1 _1 1*det 1]\.(dircos x),.y'

One other requirement is an identity matrix of appropriate length:

   id=:=@i.@#                NB. identity matrix

Now suppose the anti-clockwise angle of rotation looking outwards from the origin is t (Note that this is not according to the Right-hand-rule commonly used in mathematics, science and engineering). By a pleasing analogy with the two dimensional case the displacement components are

$-(1-\cos t)*{\vec {a}}-\sin t*{\vec {p}}$

where ${\vec {p}}$ is the cross-product between rotation axis and the point to be rotated and ${\vec {a}}$ is the result of the matrix multiplication

$${\begin{pmatrix}1-\lambda _{1}^{2}&-\lambda _{1}\lambda _{2}&-\lambda _{1}\lambda _{3}\\-\lambda _{1}\lambda _{2}&1-\lambda _{2}^{2}&-\lambda _{2}\lambda _{3}\\-\lambda _{1}\lambda _{3}&-\lambda _{2}\lambda _{3}&1-\lambda _{3}^{2}\end{pmatrix}}{\begin{pmatrix}x\\y\\z\end{pmatrix}}$$

in which the λs are the direction cosines of the axis of rotation. The parameters defining a rotation are thus an axis (three coordinates) joined to the angle t, and it seems natural to take this 4-element vector as the left argument of a rotation verb. Also many people are more comfortable with degrees rather than radians, so define:

   dtor=:180%~o.                NB. degrees to radians

rm defines the rotation matrix above and rmd multiplies it with the coordinates of the data point being rotated :

   rm =:(id - */~)@dircos        NB. rotation matrix
   rm 3 4 5
 0.82 _0.24 _0.3
_0.24  0.68 _0.4
 _0.3  _0.4  0.5
   rmd=.rm@(}:@[) +/ .* ]
   (3 4 5,dtor 60) rmd 1 2 3
_0.56 _0.08 0.4

(As an aside, taking the z axis as axis of rotation (0 0 1) so that λ₃ = 1 and λ₁ = λ₂ = 0 gives

   rm 0 0 1
1 0 0 
0 1 0 
0 0 0

which multiplies ${\textstyle {\begin{pmatrix}x\\y\\z\end{pmatrix}}}$ to give ${\textstyle {\begin{pmatrix}x\\y\\0\end{pmatrix}}}$, while the cross-product between ${\textstyle {\begin{pmatrix}0\\0\\1\end{pmatrix}}}$ and ${\begin{pmatrix}x\\y\\z\end{pmatrix}}$ is ${\textstyle {\begin{pmatrix}-y\\x\\0\end{pmatrix}}}$ so that this reduces to the formula given earlier for the two-dimensional case. )

Next define m1 and m2, bearing in mind that t has to be extracted as the 4^th element of the rotation vector :

   m1=:-.@(2&o.@({:@[))    NB. (1-cos t)
   m2=:1&o.@({:@[)    NB. sin t

Finally reflect the formula for the displacement components in

   rota=:] - (m1 * rmdata) + m2 * }:@[ xp ] NB. anti-clockwise
   rotc=:] - (m1 * rmdata) - m2 * }:@[ xp ] NB. clockwise
   (1 0 0,dtor 90) rotate 1 2 3
1 3 _2

A couple of further checks helps confirm understanding :

(a) Rotate the point (1,2,3) through a clockwise angle of 90° about the x-axis:

   (1 0 0,dtor _90)rotc 1 2 3
1 _3 2

(b) Undo a clockwise rotation of (1,2,3) with an anti-clockwise one :

   axis=.?3$10    NB. choose an axis at random
   (axis,dtor 60) rota (axis,dtor 60)rotc 1 2 3
1 2 3

Multiple data points are dealt with by, for example

   (<3 4 5,dtor 60) rota every 1 2 3;2 3 1;3 1 2
1.03505  2.5299 2.55505
3.03722 1.56268 1.52753
1.82258 0.31773 3.25227

It is noted that a similar transformation, without using displacement is described on this WikiPedia page.

Plotting star movements

Every minute of every day we all perform rotations on a merry-go-round called Earth which itself rotates continuously within an even larger solar system which itself gyrates around another even larger system and so on. This a special case of 3D rotation in which all data points in the heavens are identified by two rather than than three parameters. Astronomers measure star positions as observed from Earth in angular rather than Cartesian measure. Specifically the two angles used are altitude A which corresponds to celestial latitude, and azimuth Z which corresponds to longitude in terrestrial measurement. The stars themselves lie on the surface of a sphere called the celestial sphere which is continuously rotating about the extended Earth axis and on which every star has a latitude and a longitude which are called respectively declination d and right ascension ra. Analogous to the Greenwich meridian on Earth, the celestial sphere requires an arbitrary zero line or celestial meridian from which right ascension is measured. This is conventionally taken to be the first point in Aries, which is observable as the rightmost star in the constellation Cassiopeia. Azimuth is often measured in sidereal hours from 0 to 24 rather than degrees; the significance of ‘sidereal’ is that a sidereal year is one day longer than a solar year, that is the fixed stars appear to rotate at a slightly slower speed than the sun, the difference being about 4 minutes per day. Stars rise in the east and set in the west, and so to an Earth-bound observer looking outwards to the Pole Star, the celestial sphere appears to rotate in an anticlockwise direction.

To convert star positions defined by A and Z into (x,y,z) coordinates assume the x-axis runs west to east, the y-axis north to south, and the z-axis upward (i.e. again a left-handed coordinate system). The plane x=0 is then the celestial meridian on a fixed celestial sphere from which Z is measured clockwise (both A and Z measured in radians). (x,y,z) coordinates are then given by

$x=\cos A\sin Z;\quad y=\cos A\cos Z;\quad z=\sin A$

Inverting these formulae to convert from (x,y,z) coordinates to (A,Z) coordinates:

$A=\sin ^{-1}z;\quad Z=\cos ^{-1}{\frac {y}{\sqrt {1-z^{2}}}}(or\sin ^{-1}{\frac {x}{\sqrt {1-z^{2}}}})$ )

The diagram below shows a star S in at a point in its daily circuit around the fixed star sphere (or, as it is sometimes referred to celestial sphere):

ZSP is a spherical triangle whose sides are all great circles of the sphere and in which P is the pole and Z is the Zenith. The side SP of triangle ZSP remains constant as S proceeds along the dotted line and is equal to 90° minus the declination d. Side ZP is also constant and is equal to 90° minus the terrestrial latitude l. Quantities which change as S proceeds along its course are thus:

Side ZS = the co-latitude, that is 90° minus the altitude A

angle P = the celestial azimuth, known as the hour angle

angle Z = the terrestrial azimuth

Transits

A star is said to transit or culminate when it is at its highest point in the sky when seen by an observer on Earth. The diagram below shows a cross-section of the celestial sphere containing a star with declination 20° observed at point of transit from a latitude of 50° North which transits at (0, -cos 60°, sin 60°).

This diagram can be generalised to show that the altitude at transit is (90° - l) + d provided that d<l as in the case of the star illustrated. This star transits south, that is to the left of the zenith and dips below the terrestrial horizon for at least part of its circuit. If d>l a star is circumpolar and transits north. Here the altitude at transit is (90° + l) – d, or combining the two cases, the altitude of every star at transit is 90° – abs(l - d).

Calculating star positions

The position of a star depends on time as a parameter, which can be either local time – where was a star 6 hours ago? – or time by year – where was it 3 months ago? The star sphere appears to us to revolve from east to west completing a revolution in a sidereal day which is shorter by 1/365^th of a day (that is approximately 4 minutes) than the solar day. Thus the position of a star 6 hours ago (¼ of a day) is the same as its position 3 months ago (¼ of a year).

For example, consider the star illustrated above with declination 20^o, and ask what its position is 3 and 6 hours earlier and subsequently, that is when the hour angles is -90^o

   (0,(cos 50),(sin 50),dtor -90)rota 0,(-cos 60),sin 60
0.939693 0.219846 0.262003

This result can be confirmed by spherical trigonometry applied to the triangle ZPS. The cos formula for a spherical triangle ABC states that if a, b and c are sides measured in angles, and A, B and C are the angles between the sides with A opposite a, etc. then

${\textstyle \cos a=\cos b\cos c-\sin a\sin b\cos A}$

so applying this twice to the diagram

${\textstyle \sin d=\sin l\sin A-\cos l\cos A\cos Z}$

and

${\textstyle \sin A=\sin l\sin d-\cos l\cos d\cos H}$

Thus for the star with declination 20^o, 6 hours earlier the hour angle is -90^o so cosH = 0 and therefore

The values of a and its cosine are thus given by

    ]s=.*/sin 50 20            NB. sin a = z
0.262003
    ]cosa=.%:-.*:*/s           NB. cosa by Pythagoras
0.965067

and the azimuth value is :

   ]Z=.asin(cos 20)%cosa        NB. azimuth
76.8322

which enables the x and y coordinates to be found formulae given earlier :

   cosa*(sin,cos)Z              NB. x,y coords
0.939693 0.219846

The next step is to isolate the hour angle as a parameter (clockwise 90^o is the same as anti-clockwise -90^o):

    v=.monad :'(0,(cos 50),(sin 50),dtor y)rotc 0,(-cos 60),sin 60'

and plot values as this moves towards transit at 10^o intervals :

   v every 90 80 70 60 50 40 30 20 10 0
 0.94    0.22 0.262
0.925  0.0948 0.367
0.883 _0.0264 0.469
0.814   _0.14 0.564
 0.72  _0.243  0.65
0.604  _0.332 0.725
 0.47  _0.404 0.785
0.321  _0.457  0.83
0.163  _0.489 0.857
    0    _0.5 0.866

More generally, it is useful to convert time to angular measurement with 24 hours being equivalent to a complete rotation, which suggests a few more utility verbs :

   ttor=.o.&(%&43200@(60&#.))    NB. time (hms) to radians
   ttor 12 0 0            NB. check 12 hrs=pi rads
3.14159
   dtot=.60 60 60&#:@(*&240)    NB. deg to time (hms)
   dtot 180
12 0 0
   atod=.%&60@(60&#.)        NB. angle(deg,min)to deg
   atod 49 15
49.25

The cooordinates of the above star 15 and a half minutes after transit, are given by

   (0,(cos 50),(sin 50),ttor 0 15 30) rota 0,(-cos 60),sin 60
_0.0635044 _0.498354 0. 864645

that is a little bit to the west, a shade less south and a bit lower, all of as expected.

The next illustration concerns the sun which, unlike other stars whose declination is constant, has declination varying in the course of a year from -23.5^o to +23.5^o and back again. The sine formula in spherical trigonometry states that for a general triangle ABC :

${\textstyle {\frac {\sin a}{\sin A}}={\frac {\sin b}{\sin B}}={\frac {\sin c}{\sin C}}}$

At sunrise and sunset the sun’s altitude is zero, and so using the first of the two cosine formulae ${\textstyle \sin d=-\cos l\cos Z}$  at these times, from which ${\textstyle \cos Z=-{\frac {\sin d}{\cos l}}}$. Then using the sin formula, ${\textstyle \sin H=-{\frac {\sin X}{\cos d}}}$.

Considering London (latitude of 51^o 30') on the 21st December when the declination of the sun is -23^o 30', and the altitude at noon, that is transit, is (90^o - 51^o 30') - 23^o 30' = 15^o00' ,

   ]Z=.acos(sin 23.5)% cos atod 51 30    NB. azimuth
50.17
   ]H=.dtot asin(sin Z)%cos 23.5        NB. time to noon
3 47 27.26

Now define

cs=.(cos,sin)@atod

and use the general rotation verb rota to rotate from transit for 3 hrs 47 minutes and 27.26 seconds :

   lat=.51 30
   dec=.-23 30
   tim=.3 47 27.26
   alt=.90-atod|lat-dec
   (0,(cs lat),ttor tim)rota 0,(-cos alt),sin alt
_0.7679 _0.6405 1.278e_7

and as expected the sun is south and west at altitude zero.

Code Summary

xp=: cofs@:(norm@[ ,. ])                             NB. cross product
id=: =@i.@#                                          NB. identity matrix
norm=: % %:@(+/@:*:)                                 NB. normalise list
dircos=: % %:@(+/@:*:)                               NB. direction cosines
det=: -/ . *                                         NB. determinant
cofs=: (* signs)@:det@submats                        NB. cofactors
signs=: 1 _1&($~ #)                                  NB. successive 1 and _1
submats=: 1&(+\.)                                    NB. successive 'all but one' rows
rm=: (id - >@(*every<))@norm                         NB. rotation matrix…   id=: =@i.@#
rmdata=: rm@norm@(}:@[) +/ .* ]                      NB. …times data point
m1=: -.@(2&o.@({:@[))                                NB. multiplier (1-cos a)for rmd
m2=: 1&o.@({:@[)                                     NB. multiplier sin a for xp
rotc=: ] - (m1 * rmd) - m2 * }:@[ xp ]               NB. x=axis,angle
rota=: ] - (m1 * rmd) + m2 * }:@[ xp ]               NB. y=data point
dtor=: *&(o.%180)                                    NB. degrees to radians
dtot=: 60 60 60&#:@(*&240)                           NB. deg to time (hms)
v=: monad :'(0,(cos 50),(sin 50),dtor y)rotc 0,(-cos 60),sin 60'
ttor=: o.&(%&43200@(60&#.))                          NB. time (hms) to radians
rtod=: *&(180%(o.1))                                 NB. radians to degrees
atod=: %&60@(60&#.)                                  NB. angle(deg,min)to deg
sin=: 1&o.@dtor                                      NB. sine of an angle
cos=: 2&o.@dtor                                      NB. cosine of an angle
asin=: rtod@(_1&o.)                                  NB. arcsine of an angle
acos=: rtod@(_2&o.)                                  NB. arccosine of an angle
cs=: (cos,sin)@atod

Script

File:Fsojc48.ijs