=e==c==h=och=Calculating the intersections of two circles on a sphere.  George =Huxtable.
==h=3=20
==h=0Here are =details about my programmable-calculator program to derive the lat and =long for each crossing of two circles. It is implemented on a Casio =programmable pocket calculator, FX730P or FX795P, dating from the 1980s. =It could easily be adapted to run on other machines.

The trig routines are preset to operate in =degrees. The version of bastard-Basic that's used has just 26 possible =variables, denoted by the letters A to Z. The program lines are =numbered. Internally, it calculates in bcd (binary coded decimal) , a decimal digit at a time, rather than in true binary, and =works to a precision of 12 decimal places. 
==h=0=20
==h=0For many =functions, there is no requirement to enclose the atgument in brackets, =so for example, ACSU is the arc-cosine of variable U. A subroutine =specified by GOSUB does no more than specify an order of executing =insructions, all variables being effectively Global.
==h=0=20
==h=0Declinations =and latitudes are taken to be positive-North. Both GHAs and longitudes =are treated as increasing Westerly from Greenwich, 0 to 360 degrees, =though if preferred the final displayed long, for values > 180, could be =converted to E.
==h=0
The =program starts at line 20 with a block which reads in, via the keyboard, =the following quantities, in decimal degrees-

dec1 into B, 
==h=0GHA1 into C, =
==h=0alt1 into D, =
==h=0dec2 into E, =
==h=0GHA2 into F, =
==h=0alt2 into =G.

==h=0then as =follows, from line 80-

==h=080 =H=3DSIND*COSF/COSB-SING*COSC/COSE
90 I=3D TANE*COSC-TANB*COSF
=100 J=3DCOSF*SINC-COSC*SINF
110 =K=3DSIND*SINF/COSB-SING*SINC/COSE
120 L=3DTANE*SINC-TANB*SINF
=130 M=3D -H*I -K*L
140 N=3D =(H*I+K*L)*(H*I+K*L)-(I*I+L*L+J*J)*(H*H+K*K-J*J)		Discriminant of =quadratic in (sin=20lat).
145 IF N<0;PRINT"NO SOL": GOTO 20
150 =O=3DI*I+L*L+J*J
160 P=3D(M+SQRN)/O:IF ABSP>1;GOTO330			=	Northmost solution for sin lat1 on first pass.
165 =P=3DASNP:W=3D1:GOSUB 180						This is lat1; go =to block at 180 to find long1.
170 P=3D(M-SQRN)/O:IF ABS(P)>1;GOTO =330			 	Southmost solution for sin lat2 on second =pass.
175 P=3DASNP:W=3D-1:GOSUB180						=With lat2; to block at 180 again, for long2.
177 GOTO 20		=						All done and back to starting point.
=180 S=3DSIND/(COSP*COSB)-TANB*TANP:IF ABS(S)>1;GOTO 320	Start of =block to find long.
185 R=3D360*FRAC(1+(C+W*ACSS)/360)				=2 alternative longs for circle 1.W switches scan  
195 =S=3D360*FRAC(1+(C-W*ACSS)/360)					order on 2nd =pass.
200 U=3DSING/(COSP*COSE)-TANE*TANP:IF ABS(U)>1;GOTO 320
=205 T=3D360*FRAC(1+(F+W*ACSU)/360)					2 alternative longs for circle 2. 
215 U=3D360*FRAC(1+(F-W*ACSU)/360)
=230 IF ABS(ABS(R-T)-180)>179.9999;V=3DR:GOTO 280			These =four lines look at longs of each pair of
240 IF =ABS(ABS(R-U)-180)>179.9999;V=3DR:GOTO 280			candidate =intersections with the lat. parallel to 
250 IF ABS(ABS(S-T)-180)>179.9999;V=3DS:GOTO 280			test which =pairs are nearly the same or (just
260 IF =ABS(ABS(S-U)-180)>179.9999;V=3DS:GOTO 280			perhaps) nearly =360 different.
270 PRINT "LONG ERROR";GOTO 20
280IF P<0;GOTO =300							If Southern hemisphere.
290 PRINT RND(ABS(P),-5);"N";RND(V,-5);"W":GOTO310		If North hemi, =print lat north and westerly long.
300 PRINT =RND(ABS(P),-5);"S";RND(V,-5);"W"				If South hemi, print =lat south and westerly long.
310 RETURN
320 PRINT "LONG =ERROR2":GOTO 20
330 PRINT "LAT ERROR":GOTO 20

==h=0This =calculation follows the general lines of the method described in =sections 1 and 2 of the paper “The K-Z position solution for the =double sight”, by K H Zevering, in European Journal of Navigation, =vol 1 no 3, dec 2003. That article contained  many typographical errors which were corrected in a later issue, vol 2, no 2. =Readers should be aware that other parts of that article contain serious =errors of substance by the author, on top of the typographical errors, =but sections 1 and 2, after typos are corrected, appear to be sound. However, my treatment here uses a different =notation, to correspond with the naming of variables of my pocket =calculator.
==h=0=20
==h=0The Earth is =treated as spherical. Latitudes are North-positive. Longitudes are taken =to be positive Westwards, 0 to 360 degrees, to conform with the =direction of Hour Angles. 
==h=0=20
==h=0A celestial =body B1, at dec1, GHA1, is observed at an altitude of alt1. In which =case a circle can be drawn on a globe, centred at C1 (dec1, GHA1), with =radius z1, where zenith angle z1 =3D 90 – alt1. That circle =represents the locus of possible positions of the observer.
==h=0=20
==h=0At that same =moment, another celestial body, B2,at (dec2, GHA2) was observed at an =altitude of alt2. A second circle is drawn on the globe, centred on =(dec2, GHA2), with radius z2 (=3D 90 – alt2). That second circle =should be drawn so that it intersects the=20first. There will be two such intersections; choose one, and mark it as =P. That is the point for which we are aiming to find the position as  =(lat, long), from the intersection of those two circles. We will look =for both intersections.
==h=0=20
==h=0Join C1 and P = to the North pole N, and join C1 to P, with lines that are great =circles on the globe. We now have a spherical triangle (N to C1 to P), =with the three sides as 90-dec1, 90-lat1, and 90-alt1, and the subtended =angle at the pole (GHA1-long). In which case we can write, for that triangle, from the fundamental =cosine law-
==h=0=20
==h=0cos =(GHA1-long) =3D sin alt1 / cos lat cos dec1  - tan lat tan dec1  (eq.1)
==h=0=20
==h=0similarly, =join C2 to the North pole and to P, which creates a second spherical =triangle N to C2 to P, of which we can write-
==h=0=20
==h=0cos =(GHA2-long) =3D sin alt2 / cos lat cos dec2  - tan lat tan dec2  (eq.2)
==h=0=20
==h=0We need to =solve those two equations, finding a value for lat and long for which =they are both true, and that will represent the point P. In fact, there =will be two such solutions.
==h=0=20
==h=0The given =quantities, which define the particular problem, are dec1, GHA1, and =alt1, and dec2, GHA2, and alt2.
==h=0=20
==h=0Equations (1) =and (2) are simultaneous equations, in which lat and long are seriously =entangled together , and we have to eliminate one of those quantities. =We will choose to eliminate long. It will involve quite a lot of =algebra.
==h=0=20
==h=0Expand =cos(GHA1-long) as cosGHA1 cos long + sin GHA1 sin long, etc, and we get =
==h=0=20
==h=0cos GHA1 cos =long + sin GHA1 sin long =3D (sin alt1/cos dec1 – tan dec1 sin lat) / =cos lat  (eq.3)
==h=0cos GHA2 cos =long + sin GHA2 sin long =3D (sin alt2/cos dec2 – tan dec2 sin lat) / =cos lat  (eq.4)
==h=0=20
==h=0now multiply =eq.3 by cos GHA2 and from it, subtract eq.4 multiplied by cos GHA1, =which gives-
==h=0=20
==h=0(sin GHA1 cos =GHA2 – sin GHA2 cos GHA1)sin long  
==h=0=3D [(sin =alt1 cos GHA2 / cos dec1 – sin alt2 cos GHA1 / cos dec2) + (tan dec2 =cos GHA1 – tan dec1 cos GHA2) sin lat] / cos lat
==h=0=20
==h=0which we can =write as-
==h=0=20
==h=0sin long =3D =(H + I sin lat ) / J cos lat,    (eq. 5)
==h=0where-
==h=0H =3D  sin =alt1 cos GHA2 / cos dec1 – sin alt2 cos GHA1 / cos dec2
==h=0I  =3D  tan =dec2 cos GHA1 – tan dec1 cos GHA2    (the left-hand side is a =letter-I, not a 1)
==h=0J  =3D  sin =GHA1 cos GHA2 – sin GHA2 cos GHA1
==h=0These are =lines 80 to 100 of my calculator program.
==h=0=20
==h=0similarly, we =can multiply eq.3 by sin GHA2 and from it, subtract eq.4 multiplied by =sin GHA1, which gives-
==h=0=20
==h=0cos GHA1 sin =GHA2 – cos GHA2 sin GHA1) cos long
==h=0=3D [sin alt1 =sin GHA2 / cos dec1 – sin alt2 sin GHA1 / cos dec2 +(tan dec2 sin =GHA1 – tan dec1 sin GHA2) sin lat) / cos lat
==h=0=20
==h=0in the same =way, we can write that as-
==h=0=20
==h=0cos long =3D =(K + L sin lat) / -J cos lat     (eq. 6)
==h=0where-
==h=0K =3D  sin =alt1 sin GHA2 / cos dec1 – sin alt2 sin GHA1 / cos dec2
==h=0L =3D  tan =dec2 sin GHA1 – tan dec1 sin GHA2
==h=0These are =lines 110, 120 of my program.
==h=0=20
==h=0If we now =square both sides of eq.5, and also square both sides of eq.6, and add =them, this finally eliminates terms involving long, and results in -
==h=0=20
==h=0(sin =long)2 + (cos long)2 =3D =1 =3D [(H + I sin lat) / J cos lat ]2 + [(K + L =sin lat) / -Jcos lat]2.
==h=0or-
==h=0 =(J cos lat)2 ==3D (H + I sin lat)2 + =( K + L sin lat)2
==h=0or, =multiplying out the squares, and putting  1 - (sin =lat)2  for (cos lat)2  =
==h=0=20
==h=0=(I2 + L2 + =J2) (sin =lat)2 + (2*I*H + 2*K*L) sin lat + ( =H2 + K2  =-J2 ) =3D 0
==h=3=20
==h=0This is a =quadratic equation in (sin lat) which has the usual pair of solutions, =to be found in the traditional way.
==h=0=20
==h=0 =If we define 
==h=3M =3D - I*H - = K*L
==h=0N =3D ( I*H + =K*L )2 - ( I2 + =L2 + J2 )( =H2 + K2 - =J2 )
==h=3O =3D ( =I2 + L2 + =J2 )  (the left-hand side is a letter-O, not a =zero)
==h=0These are =lines 130 to 150
==h=0=20
==h=0then as long =as the discriminant N is not negative (which should be tested) then the =circles do, in fact, intersect, and there will be two solutions. Those =solutions are-
==h=0=20
==h=0lat1 =3D =arcsin( (M + sqrt(N)) / O), for the more Northerly intersection, and
==h=0lat2 =3D =arcsin( (M  - sqrt(N)) / O) for the more Southerly.
==h=0(lines 160, =165)
==h=0=20
==h=0Next the =longitudes of those two intersections must be found. Start with the more =Northerly, with lat =3D lat1.
==h=0=20
==h=0On each =circle, there are two possible points with that longitude, where the =latitude parallel lat1 intersects. One such point on circle 1 will have =the same longitude as one such point on circle 2, and that common =longitude is long1, the longitude of the=20intersection point that is being sought. The others are of no interest.
==h=0=20
==h=0Similarly, =there are intersections between the latitude line lat2 and the two =circles, which will give long2 of the common intersection in the same =way.
==h=0=20
==h=0Take the more =Northerly crossing point first, at lat1. 
==h=0Look first =for its intersections with circle1.
==h=0Substitute =the now-known value for lat1, into equation 1, to obtain cos =(GHA1-long1). 
==h=0Obtaining =such an angle from its cos gives two alternative results; either a =positive (Westerly) difference or the converse.
==h=0All we can =say is that either-
==h=0long1 =3D =GHA1 + arccos (sin alt1 / cos lat1 cos dec1 - tan lat1 tan dec1)  which =we can put into variable R, or
==h=0long1 =3D =GHA1 - arccos (sin alt1 / cos lat1 cos dec1 - tan lat1 tan dec1)   which =can go into S
==h=0(lines180 to =195)
==h=0=20
==h=0Look for the =intersection of the same parallel lat1 with circle 2, using equation =(2), in the same way, giving two more possible values, either-
==h=0long1 =3D =GHA2 + arccos (sin alt2 / cos lat1 cos dec2 - tan lat1 tan dec2)  which =we can put into variable T, or
==h=0long1 =3D =GHA2 - arccos (sin alt2 / cos lat1 cos dec2 - tan lat1 tan dec2)   which =can go into U
==h=0(lines 200 - =215)
==h=0It helps to =avoid confusion if those four angles are put into a common range of say =0 to 360 degrees positive.
==h=0=20
==h=0Now one =possible value, from the pair R and S, should match one from the pair T =and U. It's a simple matter to make that test for equality, for which I =have set a limit of .0001 degrees. Whichever pair are equal, that value =is taken to be the longitude sought, long1. To allow for special cases, longs which differ by 360 deg =are treated as equal.
==h=0(lines 230 to =260)
==h=0=20
==h=0The position =of intersection 1, as lat1, long1, is then displayed (line 290 or 300). =Then, the longitude of the more Southerly intersection, long2, is =determined in exactly the same way, knowing lat2, and that intersection =is displayed. The observer must be at one of those two positions, but he will have to use other evidence =to determine which. That will usually be easy.
==h=0=20
==h=0If the two =circles are centred on the same meridian, the latitudes of their two =crossings are the same. To ensure that both resulting longitudes get =flushed out, the order of testing for longitudes on the second scan is =reversed, by switching W from +1 to=20-1.
==h=0=20
==h=0I have no =idea whether the above procedure is an efficient one and have made no =attempts to improve on it. All I can say is that it does seem to give =correct results, when tested.
==h=0=20
==h=0=20
==h=0George =Huxtable, 10 June 06.
==h=0=george@huxtable.u-net.com or 1 Sandy Lane, Southmoor, Abingdon, Oxon, =OX13 5HX, UK. tel +44 1865 820222