181 lines
9.2 KiB
Text
181 lines
9.2 KiB
Text
[Haversine formula, for Rosetta Code.]
|
|
[EDSAC, initial orders 2.]
|
|
[Arrange the storage]
|
|
T49K P56F [L, 35 locations: library s/r R4 to read integer.]
|
|
[Can overwrite R9 once R9 is no longer needed.]
|
|
T46K P92F [N, 35 locations: library s/r P7 (print integer)]
|
|
T52K P128F [A, 33 locations: library subroutine T4 (arccos)]
|
|
T53K P162F [B, 44 locations: 1ibrary subroutine T1 (cos)]
|
|
T54K P206F [C, 10 locations: wrapper for library subroutine T1]
|
|
T51K P216F [G, 24 locations: inverse haversine function]
|
|
T45K P240F [H, 36 locations: haversine function]
|
|
T55K P276F [V, 11 locations: constants]
|
|
T47K P288F [M, 92 locations: main routine]
|
|
[----------------------------------------------------------------------------]
|
|
[Library subroutine R9 to read integers from tape at load time.]
|
|
[Must be placed at location 56; occupies 15 locations.]
|
|
T56K
|
|
GKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@
|
|
[----------------------------------------------------------------------------]
|
|
[Library subroutine M3. Prints header at load time and is then overwritten.]
|
|
[In this program, also sets teleprinter to figures.]
|
|
PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
|
|
*DISTANCES!IN!METRES@&WITH!RADIUS#
|
|
.. PK [after header, blank tape and PK (WWG, 1951, page 91)]
|
|
[----------------------------------------------------------------------------]
|
|
[Constants. Load at even address.]
|
|
E25K TV GK
|
|
E69K [at load time, call library s/r R9 to read 35-bit contants]
|
|
T#V [tell R9 where to store constants]
|
|
[0] 1800000000F [180 degrees, in units of 10^-7 degree]
|
|
[2] 900000000F [90 ditto]
|
|
[4] 16097821018F [for changing angle unit to radians]
|
|
[6] 13493037750F [pi/4 as EDSAC integer, i.e. times 2^34]
|
|
[8] 6371100# [Earth radius in metres, from MK-61/52 solution]
|
|
T10Z [resume normal loading]
|
|
[10] RF [17-bit 1/4]
|
|
[----------------------------------------------------------------------------]
|
|
[Main routine. Load at even address.]
|
|
E25K TM GK
|
|
[35-bit variables, general workspace]
|
|
[0] [2] [4] [6] [8]
|
|
[17-bit variables]
|
|
[10] [number of examples]
|
|
[11] [example number]
|
|
T12Z [variablss are not initialized; skip over them at load time]
|
|
[17-bit constants]
|
|
[12] PD [1]
|
|
[13] CF [ teleprinter colon (in figures mode)]
|
|
[14] @F [15] &F [16] K4096F [CR, LF, null]
|
|
[Enter with acc = 0]
|
|
[17] A8#V TD [pass Earth's radius to printer routine]
|
|
A19@ GN [finish header: print radius followed by CR, LF]
|
|
O14@ O15@
|
|
A23@ GL [call library s/r R4, 0D := number of examples]
|
|
AF T10@ [store (assume < 2^16)]
|
|
[Here acc = example number just done (0 first time)]
|
|
[27] S10@ E98@ [if all done, jump to exit]
|
|
A10@ A12@ T11@ [inc and store example number]
|
|
TD A11@ TF [0D := example number, extended to 35 bits]
|
|
A35@ GN O13@ [print example number followed by colon]
|
|
[Use library subroutine R4 to read data from tape to 0D.]
|
|
[Coordinates are lat1, lon1, lat2, lon2.]
|
|
A38@ GL AD U#@ T4D [lat1 to work{0} and 4D]
|
|
A43@ GC A4D T2#@ [0.5*cos(lat1) to work{1}]
|
|
A47@ GL AD T8#@ [lon1 to work{4}]
|
|
A51@ GL AD U4#@ T4D [lat2 to work{2} and 4D]
|
|
A56@ GC A4D T6#@ [0.5*cos(lat2) to work{3}]
|
|
A60@ GL AD S8#@ T4D [lon2 - lon1 to 4D]
|
|
A65@ GH [0.5*hav(lon2 - lon1) to 4D]
|
|
H2#@ V4D LD YF T4D [times cos(lat1) to 4D]
|
|
H6#@ V4D LD YF T8#@ [times cos(lat2) to work{4}]
|
|
A#@ S4#@ T4D [lat1 - lat2 to 4D]
|
|
A80@ GH A4D [0.5*hav(lat1 - lat2) to acc]
|
|
A8#@ T4D [add product from above, sum to 4D]
|
|
A85@ GG H4D [mult reg := angle/4 in radians]
|
|
V8#V L1F YF TD [times radius*4, to 0D for printing]
|
|
A92@ GN O14@ O15@ [print distance followed by CR, LF]
|
|
A11@ E27@ [acc := example numbr; loop back]
|
|
[Exit]
|
|
[98] O16@ [print null to flush teleprinter buffer]
|
|
ZF [stop]
|
|
[----------------------------------------------------------------------------]
|
|
[Haversine function, for Rosetta Code.]
|
|
[Input: 4D = angle x in range -360..360 deg,]
|
|
[expressed as integer multiple of 10^-7 deg.]
|
|
[Output: 4D = hav(x)/2 = (1 - cos(x))/4]
|
|
[Load at even address. Requires library subroutines R9 and T1.]
|
|
E25K TH GK
|
|
A3F T35@ [plant return link as usual]
|
|
[2] A4D E6@ [acc := x; jump if x >= 0]
|
|
T4D
|
|
[5] S4D [else x := -x]
|
|
[Here with acc = abs(x). Use symmetry about 180 deg.]
|
|
[6] S#V G10@ [jump if x < 180]
|
|
TD SD [else acc := 180 - x]
|
|
[10] A2#V [combine next 2 orders, which are shown as comments]
|
|
[A #V] [add 180: either restore x, or set x := 360 - x]
|
|
[Here 0 <= x <= 180. Use antisymmetry about 90 deg.]
|
|
[S 2 #V] [subtract 90]
|
|
E17@ [jump if x >= 90 degrees]
|
|
TD A5@ T30@ [save acc, plant S4D below]
|
|
AD G21@ [acc := x - 90, join common code]
|
|
[17] TD A2@ T30@ [save acc, plant A4D below]
|
|
SD [acc := 90 - x]
|
|
[21] A2#V [either restore x, or set x := 180 - x]
|
|
TD H4#V VD L4F YF [apply angle conversion factor]
|
|
[Here acc = x/2, now in radians. Pass it to cosine function]
|
|
T4D
|
|
A28@ GB [call cosine function; returns 4D := cos(x)/2]
|
|
[30] XF [(planted) either A4D or S4D]
|
|
RD YF A10V [divide by 2, round, add 1/4]
|
|
T4D [return haversine/2 in 4D]
|
|
[35] ZF [(planted) jump back to caller]
|
|
[----------------------------------------------------------------------------]
|
|
[Inverse haversine function.]
|
|
[Input: 4D = hav(x)/2]
|
|
[Output: 4D = x/4 in radians, where 0 <= x <= pi]
|
|
E25K TG GK
|
|
A3F T23@ [plant return link as usual]
|
|
A10V S4D LD YF [acc := 1/2 - hav(x) = cos(x)/2]
|
|
G13@ [jumo if cos(x) < 0]
|
|
T4D [pass cos(x)/2 to library subroutine T4 (arccis)]
|
|
A8@ GA AD RD [call T4, 0D ;= x/2; then acc := x/4]
|
|
E21@ [join common code]
|
|
[13] TD SD T4D [here if cos(x) < 0; set 4D := -cos(x)/2]
|
|
A16@ GA [call T4 (arccos); 0D := (pi - x)/2]
|
|
SD RD A6#V [acc := -(pi - x)/4 + pi/4 = x/4]
|
|
[21] YF T4D [round redult and return it in 4D]
|
|
[23] ZF [(planted) jump back to caller]
|
|
[----------------------------------------------------------------------------]
|
|
[Wrapper for library subroutine T1 (cosine).]
|
|
[Input: 4D = angle x as integer multiple of 10^-7 degree.]
|
|
[Output: 4D = cos(x)/2]
|
|
E25K TC GK
|
|
A3F T9@ [plant return link as usual]
|
|
H4#V V4D L4F YF T4D [4D := x/2 in radians]
|
|
A7@ GB [call T1, puts cos(x)/2 into 4D]
|
|
[9] ZF [jump back to caller]
|
|
[----------------------------------------------------------------------------]
|
|
[Library subroutine T1: calculates cos(x), where abs(x) <= pi/2]
|
|
[Input: 4D = x/2]
|
|
[Output: 4D = cos(x)/2]
|
|
[Requires library subroutine R9.]
|
|
E25K TB
|
|
GKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@E13ZT32#@
|
|
1614F73454F243967F54539267F763549741F5726623061#
|
|
TZA3FT30@H4DV4DYFT4DH4DN32#@A34#@TDNDA36#@TDNDA38#@
|
|
TDNDA40#@TDNDA42#@TDNDYFTDNDS4DA31@YFT4DEFIFT44Z
|
|
[----------------------------------------------------------------------------]
|
|
[Library subroutine T4: calculates arccos]
|
|
[Input: 4D = cos(x)/2, assumed >= 0]
|
|
[Output: 0D = x/2 where 0 <= x <= pi/2]
|
|
E25K TA
|
|
GKA3FT28@TDA32@T6DH4DV4DL1FS29@YFE16@T4DSDA6DTDS4DT4DA6DRD
|
|
YFG4@HDN30#@YFE27@T4DS4DTDEFIFT30#ZD888FT30ZO699DT32ZK4096F
|
|
[----------------------------------------------------------------------------]
|
|
[Library subroutine P7. Prints 35-bit positive integer in 0D.]
|
|
[Up to 10 digits, right-justified, padded on left with spaces.]
|
|
E25K TN
|
|
GKA3FT26@H28#@NDYFLDT4DS27@TFH8@S8@T1FV4DAFG31@SFLDUFOFFFSF
|
|
L4FT4DA1FA27@G11@T28#ZPFT27ZP1024FP610D@524D!FO30@SFL8FE22@
|
|
[----------------------------------------------------------------------------]
|
|
[Library subroutine R4. Reads signed integer from tape at run time.]
|
|
[Input : None]
|
|
[Output : 0D = 35-bit signed integer.]
|
|
E25K TL
|
|
GKA3FT21@T4DH6@E11@P5DJFT6FVDL4FA4DTDI4FA4FS5@G7@S5@G20@SDTDT6FEF
|
|
[----------------------------------------------------------------------------]
|
|
E25K TM GK [M parameter again]
|
|
E17Z [define entry point]
|
|
PF [acc = 0 on entry]
|
|
[----------------------------------------------------------------------------]
|
|
[Latitude and longitude of each airport, in units of 10^-7 degree.]
|
|
[Read from tape by library subroutine R4. Note that sign comes after value.]
|
|
3+ [number of examples]
|
|
361200000+866700000- [Nashville]
|
|
339400000+1184000000- [Los Angeles]
|
|
547166667+205000000+ [Kaliningrad]
|
|
431166667+1319000000+ [Vladivostok]
|
|
499133333+62916667- [St Mary's, Scilly Isles]
|
|
601919444+12436111- [Tingwall, Shetland Isles]
|