A reconstruction of VIEW, the Manned Spacecraft Center program that drew Apollo window views on a UNIVAC 1108 (FORTRAN V, output photographed from a CRT onto microfilm: NASA TN D-6853, p. 3). We know of no surviving source code; the geometry here is our own FORTRAN and this page only films it.
src/vdrive.f
C=======================================================================
C
C V I E W - 1 1 0 8 WINDOW VIEW KERNEL: FRAME DRIVER
C
C Draws what an Apollo 11 crewman would see out of a window, as
C line vectors and star points in plot degrees, for a film
C recorder to expose. After the MSC program VIEW (G. B. Roush;
C documented by A. N. Lunde and C. T. Hyle). New code; the target
C is the surviving output, MSC IN 69-FM-197 and the film clip.
C
C ENTRY POINTS (called by the chassis, shell.f90; also SIMRUN in
C sim.f, which runs the engine and fills the tape)
C VINIT (ISC, GET, YAW, PIT, ROL, FOV)
C select scene ISC and return its default inputs.
C VFRAME (GET, YAW, PIT, ROL, FOV, IFLAG,
C VB, NV, SB, NS, LB, NL, HD, TB, NT, TC, NCH)
C draw one frame. VB(5,MAXV) line vectors X1 Y1 X2 Y2
C STYLE, SB(3,MAXS) points X Y MAG, LB(4,MAXL) labels
C X Y KIND ID, HD(24) header values, TB/TC text records.
C
C ELEMENTS. The kernel is a set of separately compiled elements,
C linked by the build, as an 1108 program was put together by
C the Collector, which "is a system processor designed to provide
C the user with a means of gathering (collecting) and
C interconnecting one or more relocatable elements to produce a
C program" (UE-637 sec. 5.1; docs/batch-pipeline.md).
C vdrive.f this driver: scenes, cameras, model placement
C vlayer.f the layer dispatcher and each scene's layer list
C Core: ephem.f (time, Sun, Moon), traj.f (trajectory legs,
C the replay), sim.f (the engine), tape.f (the tape it
C writes), vsrc.f (the state source: replay or tape),
C vview.f (camera target and external view),
C pen.f (projection, clipping, visibility, vectors),
C vtext.f (text records), vmath.f (vectors, matrices),
C models.f (spacecraft model library)
C Layers, one per drawable, all called as
C LAYER(GET, VB, NV, SB, NS, LB, NL) by LAYERS:
C lframe.f 1 plot frame, lstars.f 2 stars, lsun.f 3 Sun,
C lmoon.f 4 Moon and craters (lmoon6.f its whole-disc
C extras), learth.f 5 Earth, lvehic.f 6 vehicles
C (lvlab.f their labels and markers),
C lcoas.f 7 COAS reticle, lshad.f 8 LM shadow,
C llpd.f 9 LPD and LM window
C Data: viewdata.f (BLOCK DATA, generated), viewcom.inc COMMON
C
C THE ELEMENTS AGAINST TN D-6853, printed p. 3 (our reading). "The
C program consists of two basic parts: the integrator portion and
C the graphic-display portion"; the Apollo modifications "were
C associated with the input/output options, coordinate
C transformations, lunar- and solar-ephemeris installation, three-
C dimensional-display problems, and realistic spacecraft-window
C outlines".
C integrator portion sim.f (and traj.f's replay)
C graphic-display portion pen.f, the layers via vlayer.f
C ephemeris installation ephem.f
C coordinate transformations vmath.f, the frames in vdrive.f
C three-dimensional display pen.f, models.f
C window outlines window and cabin models (to come)
C input/output vtext.f, the plot-tape buffers,
C the scenarios (data/scenarios),
C the tape (tape.f)
C OUR READING: THE REPORT NAMES FUNCTIONS, NOT FILES.
C
C PROJECTION
C Radially symmetric about the boresight, gnomonic to
C stereographic with the field; see PROJ (pen.f).
C
C=======================================================================
SUBROUTINE VINIT(ISC, GET, YAW, PIT, ROL, FOV)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER ISC
DOUBLE PRECISION GET, YAW, PIT, ROL, FOV
DOUBLE PRECISION R(3), V(3), PM(3), E(3), S(3), X, Y, TFIX
INTEGER I, ISNSC(9)
DOUBLE PRECISION VDOT, EVGET
C The scenario of each scene: Apollo 11 as flown for scenes 1-8,
C Apollo 8 as flown for scene 9.
DATA ISNSC / 1, 1, 1, 1, 1, 1, 1, 1, 2 /
C
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (INITD .NE. 1) THEN
CALL TABSET
CALL MLIB
ISN = 0
INITD = 1
END IF
C RESTOMOD END
ISCN = ISC
IF (ISCN .LT. 1 .OR. ISCN .GT. 9) ISCN = 1
IF (ISNSC(ISCN) .NE. ISN) CALL SNSET(ISNSC(ISCN))
YAW = 0.0D0
PIT = 0.0D0
ROL = 0.0D0
C
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (ISCN .EQ. 1 .OR. ISCN .EQ. 9) THEN
C EARTHRISE. Scene 1: one minute before the Earth's disc
C clears the lunar horizon on the revolution before the landing.
C The view is turned in azimuth (AZOFF) so the Earth rises
C mid-frame.
C Scene 9, APOLLO 8 EARTHRISE: at the photograph AS08-14-2383
C (the scenario's PHOTO event), with the field of its 250 mm lens
C on the 70 mm frame, 12.72 deg: the lens from the Apollo 8
C Flight Journal (day 4, orbit 4, commentary: "another Hasselblad
C with a 250-mm lens"), the 55.74 mm gate measured by us on the
C ASU scan of AS08-14-2383 (tothemoon.im-ldi.com), taking the
C 70 mm perforation pitch as 4.75 mm (unsourced). The azimuth is
C set at the photograph's time; the pointing is S9REF's.
CALL ERFIND
GET = TERISE - 60.0D0
FOV = 8.0D0
IF (ISCN .EQ. 9) GET = EVGET(KEPHO)
IF (ISCN .EQ. 9) FOV = 12.72D0
CALL VSTATE(GET, 2, R, V)
CALL MOONG(GET, PM)
E(1) = -PM(1) - R(1)
E(2) = -PM(2) - R(2)
E(3) = -PM(3) - R(3)
CALL VUNIT(R)
CALL VUNIT(V)
CALL VCRS(V, R, S)
X = VDOT(E, V)
Y = VDOT(E, S)
AZOFF = DATAN2(Y, X) / DR
ELSE IF (ISCN .EQ. 2) THEN
C EARTH APPROACH on the transearth coast. The attitude is
C held inertially: boresight ELOFF deg ahead of the Earth's
C centre as seen at TFIX, up against the Earth's drift across
C the sky, so the disc climbs in from below, grows, and leaves
C only its limb arc as the spacecraft closes on entry.
FXDT = 10.0D0 * 3600.0D0
TFIX = TETP - FXDT
GET = TETP - 5.0D0 * 3600.0D0
FOV = 60.0D0
ELOFF = 0.0D0
CALL VSTATE(TFIX, 1, R, V)
CALL VSTATE(TETP - 3600.0D0, 1, E, S)
DO 22 I = 1, 3
PM(I) = -R(I)
E(I) = -E(I)
22 CONTINUE
CALL VUNIT(PM)
CALL VUNIT(E)
C Net drift of the Earth centre across the line of sight from
C TFIX to an hour before entry.
X = VDOT(E, PM)
DO 24 I = 1, 3
S(I) = E(I) - X * PM(I)
24 CONTINUE
CALL VUNIT(S)
Y = ELOFF * DR
DO 26 I = 1, 3
FXB(I) = DCOS(Y) * PM(I) + DSIN(Y) * S(I)
FXU(I) = -(DCOS(Y) * S(I) - DSIN(Y) * PM(I))
26 CONTINUE
ELSE IF (ISCN .EQ. 3) THEN
C EARTH PARKING ORBIT, looking forward at the horizon.
GET = 1.5D0 * 3600.0D0
FOV = 70.0D0
ELSE IF (ISCN .EQ. 4) THEN
C LM RENDEZVOUS / INSPECTION, two minutes after undocking.
GET = EVGET(KEUND) + 120.0D0
FOV = 12.0D0
ELSE IF (ISCN .EQ. 8) THEN
C THE DOCKED STACK IN TRANSLUNAR COAST (a modern addition: VIEW
C drew vehicles as seen from a vehicle, TN D-6853 p. 12, not
C from outside both). Half an hour into passive thermal
C control, which began at 10:58:19 (the scenario's PTC event),
C after the LM's extraction (4:17) and before the first
C midcourse correction (26:45): our choice of moment.
GET = EVGET(KEPTC) + 1800.0D0
FOV = 40.0D0
ELSE IF (ISCN .EQ. 7) THEN
C TRANSPOSITION AND DOCKING. Mid-approach, about 56 ft out
C (see S7POSE for the closing law), looking along the CSM +X
C axis at the LM docking target. Field of view: ours.
GET = EVGET(KEAPR) + 60.0D0
FOV = 30.0D0
ELSE IF (ISCN .EQ. 6) THEN
C MOON VIEW (a modern addition, not a 1969 plot type we have a
C source for). The camera sits 35,000 km above the sub-observer
C point, our choice, which with the field below puts the disc at
C 85 percent of the frame; yaw and pitch then move the
C sub-observer point (see SCNCAM). GET at touchdown, so the
C terminator falls as it did for the landing.
GET = LUT0
S6DST = RM + 35000.0D0
FOV = DBLE(NINT(20.0D0 * DASIN(RM / S6DST) / DR / 0.85D0))
& / 10.0D0
ELSE
C LM DESCENT, commander's front window, P64 approach.
GET = 102.0D0*3600.0D0 + 42.0D0*60.0D0
C The film's descent frame is numbered to +-50 at its edges; in
C the gnomonic plot (see PROJ) that is a physical field of
C 2 ATAN(50 deg in radians) = 82.4 deg.
FOV = 82.4D0
END IF
C RESTOMOD END
RETURN
END
C
C=======================================================================
SUBROUTINE VFRAME(GET, YAW, PIT, ROL, FOV, IFLAG,
& VB, NV, SB, NS, LB, NL, HD, TB, NT, TC, NCH)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, YAW, PIT, ROL, FOV
INTEGER IFLAG, NV, NS, NL
DOUBLE PRECISION VB(5,MAXV), SB(3,MAXS), LB(4,MAXL), HD(24)
DOUBLE PRECISION TB(4,MAXT)
INTEGER NT, TC(MAXTC), NCH
DOUBLE PRECISION PM(3), CG(3), CV(3), RB, RNG, D1, D2, D3, D4
DOUBLE PRECISION PB(3), RR, VNRM, VDOT, RHO
INTEGER I, IREF, IWIN, IOK, J, LOOKD
DOUBLE PRECISION MR1(3,3), MR2(3,3)
C
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (INITD .NE. 1 .OR. ISCN .EQ. 0) THEN
CALL VINIT(1, D1, D2, D3, D4, RB)
END IF
C RESTOMOD END
NV = 0
NS = 0
NL = 0
DO 10 I = 1, 24
HD(I) = 0.0D0
10 CONTINUE
IFLG = IFLAG
C Label level (VSETIN). With in_lablv 0, in_flags bit 0 means all
C labels, as before; with 1 or more the level decides and bit 0
C is set here, so the names are lettered (TXALL).
ILEV = 3 * MOD(IFLG, 2)
IF (ILABL .GE. 1) ILEV = ILABL
IF (ILABL .GE. 1 .AND. MOD(IFLG, 2) .EQ. 0) IFLG = IFLG + 1
C State source: in_flags bit 3, the tape if the engine has run.
ISRC = MOD(IFLG / 8, 2)
ISRCU = 0
FOVH = 0.5D0 * FOV
IF (FOVH .LT. 0.05D0) FOVH = 0.05D0
IF (FOVH .GT. 89.0D0) FOVH = 89.0D0
C Projection constant (see PROJ), box half-width in plot deg, the
C angle to the frame corner, and the projection's usable limit.
PK = 1.0D0 + DMIN1(1.0D0, DMAX1(0.0D0, (2.0D0*FOVH - 100.0D0)
& / 70.0D0))
BOXH = PK * DTAN(FOVH * DR / PK) / DR
THVIEW = PK * DATAN(1.4143D0 * BOXH * DR / PK) / DR + 0.5D0
THLIM = 90.0D0 * PK - 0.5D0
IF (THVIEW .GT. THLIM) THVIEW = THLIM
CSVIEW = DCOS(THVIEW * DR)
C
C World at this GET.
CALL TSET(GET)
CALL MOONG(GET, PM)
CALL SUNG(GET, SUNU)
CALL MOONRT(GET, MMF)
C Earth fixed to J2000: GMST about the pole of date, then the
C precession back to J2000 (PRECM).
CALL ROTZ(GMST, MR1)
CALL PRECM(TCEN, MR2)
DO 18 I = 1, 3
DO 17 J = 1, 3
MEF(I,J) = MR2(1,I) * MR1(1,J) + MR2(2,I) * MR1(2,J)
& + MR2(3,I) * MR1(3,J)
17 CONTINUE
18 CONTINUE
C
C Camera position CG (geocentric), velocity CV relative to the
C reference body IREF, window code IWIN, reference attitude.
S6LAT = PIT
S6LON = YAW
CALL SCNCAM(GET, PM, CG, CV, IREF, IWIN)
C Spacecraft models first: they hide stars and bodies.
CALL SCNMOD(GET)
C The camera target and the external view (vview.f).
CALL VIEWPT(GET, PM, CG, YAW, PIT, ROL, LOOKD)
DO 20 I = 1, 3
EPOS(I) = -CG(I)
MPOS(I) = PM(I) - CG(I)
20 CONTINUE
CALL MTXV(MMF, MPOS, CAMF)
DO 30 I = 1, 3
CAMF(I) = -CAMF(I)
30 CONTINUE
C
C Free look, then the boresight in the Moon frame.
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (ISCN .EQ. 6) THEN
CALL LOOK(0.0D0, 0.0D0, ROL)
ELSE IF (LOOKD .EQ. 0) THEN
CALL LOOK(YAW, PIT, ROL)
END IF
C RESTOMOD END
CALL MTXV(MMF, CB, CBMF)
C
ISTYLE = 1
C The scene's layers, in order (vlayer.f).
CALL LAYERS(GET, VB, NV, SB, NS, LB, NL)
C
C Header.
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IREF .EQ. 1) THEN
RNG = VNRM(EPOS)
RB = RE
ELSE
RNG = VNRM(MPOS)
RB = RM
END IF
C RESTOMOD END
HD(1) = GET
HD(2) = 2.0D0 * FOVH
HD(3) = RNG / 1.852D0
HD(4) = (RNG - RB) / 1.609344D0
HD(5) = VNRM(CV) * 1000.0D0 / 0.3048D0
HD(6) = DBLE(IREF)
HD(7) = DBLE(ISCN)
HD(8) = DBLE(IWIN)
IF (ISCN .EQ. 4) HD(9) = 300.0D0
C Scene 5: the footpads' altitude (LMDESC), 0 at touchdown.
IF (ISCN .EQ. 5) HD(10) = LMALT * 1000.0D0 / 0.3048D0
IF (ISCN .EQ. 7) HD(9) = S7RNG
C Reference body in the picture, for the page's camera steering:
C centre X, Y (deg, even off frame), angular radius, in front flag.
C Scene 4: the LM. Scene 7: the LM's docking target (S7POSE).
DO 40 I = 1, 3
PB(I) = EPOS(I)
IF (IREF .EQ. 2) PB(I) = MPOS(I)
IF (ISCN .EQ. 4) PB(I) = 300.0D0 * 0.3048D-3 * BREF(I)
IF (ISCN .EQ. 7) PB(I) = S7LP(I) - 0.72D-3 * S7AT(I,2)
IF (ISCN .EQ. 8) PB(I) = MDP(I,KCSM) + 3.2D-3 * MDAT(I,1,KCSM)
40 CONTINUE
RR = RE
IF (IREF .EQ. 2) RR = RM
IF (ISCN .EQ. 4) RR = 4.5D-3
C Scene 7: the LM, half its 14 ft 1 in width (Apollo 11 press kit,
C printed p. 96).
IF (ISCN .EQ. 7) RR = 2.15D-3
C Scene 8: about half the stack's length.
IF (ISCN .EQ. 8) RR = 1.0D-2
CALL PROJ(PB, HD(11), HD(12), IOK)
HD(13) = RHO(DASIN(DMIN1(1.0D0, RR / VNRM(PB))))
HD(15) = BOXH
C The scenario's epoch as an offset (s) from Apollo 11 range zero,
C for the page's clock: UTC = 1969-07-16 13:32:00 + HD(16) + GET.
HD(16) = (TJD0 - JD0) * 86400.0D0
C The state source used (0 replay, 1 sim corrected, 2 sim free),
C and the last engine run's error at the reference row nearest
C GET: position (km), velocity (ft/s), that row's g.e.t. (s).
HD(17) = DBLE(ISRCU)
CALL SIMERR(GET, J, HD(20), HD(18), HD(19))
HD(14) = 0.0D0
IF (VDOT(PB, CB) .GT. 0.0D0) HD(14) = 1.0D0
C The vehicles in this frame's world (VPRES): 1 CSM, 2 LM, 4 S-IVB.
CALL VPRES(GET)
HD(21) = DBLE(IVBIT)
C Text for the recorder's character generator.
CALL TXALL(LB, NL, TB, NT, TC, NCH)
RETURN
END
C
C=======================================================================
C SCENE CAMERAS. Position, velocity, reference body and the
C reference attitude RREF, UREF, BREF for scene ISCN at GET.
C=======================================================================
SUBROUTINE SCNCAM(GET, PM, CG, CV, IREF, IWIN)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, PM(3), CG(3), CV(3)
INTEGER IREF, IWIN
DOUBLE PRECISION R(3), V(3), RU(3), VU(3), SU(3), H(3), E(3)
DOUBLE PRECISION DIP, CA, SA, CD, SD, P2(3), R2(3), K, VDOT
DOUBLE PRECISION XB(3), YB(3), ZB(3), PMF(3), VNRM
INTEGER I
C
IWIN = 1
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (ISCN .EQ. 1 .OR. ISCN .EQ. 4 .OR. ISCN .EQ. 9) THEN
C CSM in lunar orbit.
IREF = 2
CALL VSTATE(GET, 2, R, V)
DO 10 I = 1, 3
CG(I) = PM(I) + R(I)
CV(I) = V(I)
RU(I) = R(I)
VU(I) = V(I)
10 CONTINUE
CALL VUNIT(RU)
CALL VUNIT(VU)
IF (ISCN .EQ. 1 .OR. ISCN .EQ. 9) THEN
C Forward along the orbit, turned AZOFF in azimuth, down to
C the horizon by the dip angle.
CALL VCRS(VU, RU, SU)
CA = DCOS(AZOFF * DR)
SA = DSIN(AZOFF * DR)
DIP = DACOS(RM / VNRM(R))
CD = DCOS(DIP)
SD = DSIN(DIP)
DO 20 I = 1, 3
H(I) = CA * VU(I) + SA * SU(I)
BREF(I) = CD * H(I) - SD * RU(I)
UREF(I) = SD * H(I) + CD * RU(I)
20 CONTINUE
C Scene 9: turned to the photograph's framing (S9REF).
IF (ISCN .EQ. 9) CALL S9REF
ELSE
C Out of plane toward the LM, local vertical up.
CALL VCRS(RU, VU, BREF)
DO 30 I = 1, 3
UREF(I) = RU(I)
30 CONTINUE
END IF
ELSE IF (ISCN .EQ. 2) THEN
C Coast. Attitude held inertially (FXB, FXU, set by VINIT).
IREF = 1
CALL VSTATE(GET, 1, R, V)
DO 40 I = 1, 3
CG(I) = R(I)
CV(I) = V(I)
BREF(I) = FXB(I)
UREF(I) = FXU(I)
40 CONTINUE
ELSE IF (ISCN .EQ. 3) THEN
C Parking orbit. Forward, boresight 8 deg above the horizon.
IREF = 1
CALL VSTATE(GET, 1, R, V)
DO 70 I = 1, 3
CG(I) = R(I)
CV(I) = V(I)
RU(I) = R(I)
VU(I) = V(I)
70 CONTINUE
CALL VUNIT(RU)
CALL VUNIT(VU)
DIP = DACOS(RE / VNRM(R)) - 8.0D0 * DR
CD = DCOS(DIP)
SD = DSIN(DIP)
DO 80 I = 1, 3
BREF(I) = CD * VU(I) - SD * RU(I)
UREF(I) = SD * VU(I) + CD * RU(I)
80 CONTINUE
ELSE IF (ISCN .EQ. 6) THEN
C Moon view: looking at the centre from above the sub-observer
C point (S6LAT, S6LON), selenographic north up (our choice; a
C J2000-north layout would do as well).
IREF = 2
E(1) = DCOS(S6LAT * DR) * DCOS(S6LON * DR)
E(2) = DCOS(S6LAT * DR) * DSIN(S6LON * DR)
E(3) = DSIN(S6LAT * DR)
CALL MXV(MMF, E, H)
DO 95 I = 1, 3
CG(I) = PM(I) + S6DST * H(I)
CV(I) = 0.0D0
BREF(I) = -H(I)
UREF(I) = MMF(I,3)
95 CONTINUE
K = VDOT(UREF, BREF)
DO 96 I = 1, 3
UREF(I) = UREF(I) - K * BREF(I)
96 CONTINUE
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (VDOT(UREF, UREF) .LT. 1.0D-12) THEN
DO 97 I = 1, 3
UREF(I) = MMF(I,1)
97 CONTINUE
END IF
C RESTOMOD END
CALL VUNIT(UREF)
ELSE IF (ISCN .EQ. 8) THEN
C The docked stack (S8ATT) from 60 m on its far side from the
C Earth, looking at it with the Earth behind; the stack's X axis
C up. Its own view is external (VIEWPT), which yaw and pitch
C carry around the stack.
IREF = 1
CALL VSTATE(GET, 1, R, V)
CALL S8ATT(GET)
DO 120 I = 1, 3
BREF(I) = -R(I)
120 CONTINUE
CALL VUNIT(BREF)
K = VDOT(S8AT(1,1), BREF)
DO 125 I = 1, 3
UREF(I) = S8AT(I,1) - K * BREF(I)
CG(I) = R(I) - 0.060D0 * BREF(I)
CV(I) = V(I)
125 CONTINUE
CALL VUNIT(UREF)
ELSE IF (ISCN .EQ. 7) THEN
C Transposition and docking. The CSM on the translunar ellipse
C (30 m from the S-IVB, nothing at this scale), turned around to
C face the stack, which holds an inertial attitude (S7ATT).
C Boresight along the CSM +X axis, which is down the LM's -X
C axis; the LM front (+Z) up.
IREF = 1
CALL VSTATE(GET, 1, R, V)
CALL S7ATT
DO 110 I = 1, 3
CG(I) = R(I)
CV(I) = V(I)
BREF(I) = -S7AT(I,1)
UREF(I) = S7AT(I,3)
110 CONTINUE
ELSE
C LM descent. Boresight 46 deg down from the LM +Z axis in
C the X-Z plane, so the LPD scale 0..80 deg runs top to bottom.
IREF = 2
IWIN = 2
CALL LMDESC(GET, PMF, XB, YB, ZB)
CALL LMDESC(GET + 1.0D0, P2, H, E, R2)
CALL MXV(MMF, PMF, R)
CALL MXV(MMF, P2, R2)
DO 90 I = 1, 3
CG(I) = PM(I) + R(I)
CV(I) = R2(I) - R(I)
90 CONTINUE
CALL MXV(MMF, XB, H)
CALL MXV(MMF, ZB, E)
CD = DCOS(46.0D0 * DR)
SD = DSIN(46.0D0 * DR)
DO 100 I = 1, 3
BREF(I) = CD * E(I) - SD * H(I)
UREF(I) = SD * E(I) + CD * H(I)
100 CONTINUE
END IF
C RESTOMOD END
CALL VCRS(BREF, UREF, RREF)
CALL VUNIT(RREF)
RETURN
END
C
C-----------------------------------------------------------------------
C FREE LOOK. Camera axes from the reference attitude turned by
C yaw (+ right), pitch (+ up) and roll (+ camera clockwise as the
C viewer sees it, so the picture turns counter-clockwise).
C-----------------------------------------------------------------------
SUBROUTINE LOOK(YAW, PIT, ROL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION YAW, PIT, ROL
DOUBLE PRECISION B1(3), R1(3), U2(3), C, S
INTEGER I
C = DCOS(YAW * DR)
S = DSIN(YAW * DR)
DO 10 I = 1, 3
B1(I) = C * BREF(I) + S * RREF(I)
R1(I) = C * RREF(I) - S * BREF(I)
10 CONTINUE
C = DCOS(PIT * DR)
S = DSIN(PIT * DR)
DO 20 I = 1, 3
CB(I) = C * B1(I) + S * UREF(I)
U2(I) = C * UREF(I) - S * B1(I)
20 CONTINUE
C = DCOS(ROL * DR)
S = DSIN(ROL * DR)
DO 30 I = 1, 3
CU(I) = C * U2(I) + S * R1(I)
CR(I) = C * R1(I) - S * U2(I)
30 CONTINUE
RETURN
END
C
C=======================================================================
C TABLES. Crater centres and coastline points to unit vectors.
C=======================================================================
SUBROUTINE TABSET
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION A, B
INTEGER I
DO 10 I = 1, NCRAT
A = CRLAT(I) * DR
B = CRLON(I) * DR
CRV(1,I) = DCOS(A) * DCOS(B)
CRV(2,I) = DCOS(A) * DSIN(B)
CRV(3,I) = DSIN(A)
10 CONTINUE
DO 20 I = 1, NCPT
A = CLAT(I) * DR
B = CLON(I) * DR
CEV(1,I) = DCOS(A) * DCOS(B)
CEV(2,I) = DCOS(A) * DSIN(B)
CEV(3,I) = DSIN(A)
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SCNMOD: place this frame's models (after SCNCAM, before the sky
C is drawn, since placed solids hide stars and bodies).
C-----------------------------------------------------------------------
SUBROUTINE SCNMOD(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET
CALL MCLEAR
IF (ISCN .EQ. 4) CALL LMPIRO(GET)
IF (ISCN .EQ. 7) CALL S7POSE(GET)
IF (ISCN .EQ. 8) CALL S8POSE
RETURN
END
C
C-----------------------------------------------------------------------
C LMPIRO: the LM 300 ft from the CSM along the reference
C boresight, turning slowly for inspection (scene 4).
C-----------------------------------------------------------------------
SUBROUTINE LMPIRO(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET
DOUBLE PRECISION AT(3,3), R1(3,3), R2(3,3), R3(3,3), R4(3,3)
DOUBLE PRECISION BX(3,3), LP(3), BO(3), T, PS, TH, PH, DIST
DOUBLE PRECISION EVGET
INTEGER I
C Body axes in the reference frame: X up, Z toward the camera.
DO 10 I = 1, 3
BX(I,1) = UREF(I)
BX(I,2) = -RREF(I)
BX(I,3) = -BREF(I)
10 CONTINUE
T = GET - EVGET(KEUND)
PS = (35.0D0 + 1.5D0 * T) * DR
TH = -(20.0D0 + 15.0D0 * DSIN(T * 2.0D0 * PI / 300.0D0)) * DR
PH = (10.0D0 * DSIN(T * 2.0D0 * PI / 420.0D0)) * DR
C Yaw about body X, tilt toward the camera about body Y, then
C turn in the picture about body Z (toward the camera).
CALL ROTX(PS, R1)
CALL ROTY(TH, R2)
CALL ROTZ(PH, R3)
CALL MXM(R3, R2, R4)
CALL MXM(R4, R1, R2)
CALL MXM(BX, R2, AT)
DIST = 300.0D0 * 0.3048D-3
DO 20 I = 1, 3
LP(I) = DIST * BREF(I)
20 CONTINUE
C Centred on the stage joint.
CALL SETV(BO, 2.3D0, 0.0D0, 0.0D0)
CALL MPLACE(KLMD, AT, LP, BO)
RETURN
END
C
C=======================================================================
C TRANSPOSITION AND DOCKING (scene 7). TN D-6853 (printed p. 12)
C lists among VIEW's capabilities integrating trajectories "after
C separation of the LM and the CSM or the CSM/LM and S-IVB
C vehicles", drawing "the vehicle outlines of the CSM, LM, and the
C S-IVB" at their apparent size, and "hidden-line models of the LM
C and the S-IVB". The contents and OCR of the Apollo 11 report
C (MSC IN 69-FM-197) list no such views: no answer key.
C
C Apollo 11 times (Mission Report MSC-00171, table 7-II, printed
C p. 7-9): command module/S-IVB separation 3:17:04.6, docking
C 3:24:03.1. The crew expected to be "out about 66 feet", and
C Collins guessed "around 100 or so" (Apollo 11 Flight Journal,
C 003:38:07); the report has "a maximum separation distance of at
C least 100 feet" and contact "at an estimated 0.1 ft/sec"
C (printed p. 4-2). Our closing law: 100 ft until 3:20:30 (the
C turnaround is not modelled; the time is our guess, between the
C separation and Collins' "I'm still quite a ways" at 3:22:25),
C then range = 100 (A U + (1-A) U**2) ft, U the fraction of the
C closing time left, A set so the range rate at contact is 0.1
C ft/s.
C-----------------------------------------------------------------------
SUBROUTINE S7POSE(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, TCLS, TDOK, U, A, D, V(3), W(3), Z(3)
DOUBLE PRECISION PS(3), EVGET
INTEGER I
C Approach start and docking from the scenario (APPR, DOCK).
TCLS = EVGET(KEAPR)
TDOK = EVGET(KEDOK)
U = (TDOK - GET) / (TDOK - TCLS)
IF (U .GT. 1.0D0) U = 1.0D0
IF (U .LT. 0.0D0) U = 0.0D0
A = 0.1D0 * (TDOK - TCLS) / 100.0D0
S7RNG = 100.0D0 * (A * U + (1.0D0 - A) * U * U)
C The camera is the COAS in the CSM's left rendezvous window,
C looking parallel to the docking axis 0.72 m to the LM's -Y side
C of it and 2 m behind the CSM's docking ring (our guesses; the
C target is set off to match, LMDOCK). S7BO: the tunnel top.
CALL SETV(S7BO, 4.35D0, 0.0D0, -0.6D0)
D = (S7RNG * 0.3048D0 + 2.0D0) * 1.0D-3
DO 10 I = 1, 3
S7LP(I) = D * BREF(I) + 0.72D-3 * S7AT(I,2)
10 CONTINUE
CALL MPLACE(KLMS, S7AT, S7LP, S7BO)
C The S-IVB on the same axes, the top of its IU 1.5 m below the
C LM's base on the descent stage's axis (Y = Z = 0), which puts
C the LM tunnel near the top of the 28 ft SLA (press kit, printed
C p. 88) with room for the SPS nozzle: our guess. Our LM model
C has its tunnel 0.6 m aft of that axis.
CALL SETV(V, (-1.5D0 - S7BO(1)) * 1.0D-3, -S7BO(2) * 1.0D-3,
& -S7BO(3) * 1.0D-3)
CALL MXV(S7AT, V, W)
DO 20 I = 1, 3
PS(I) = S7LP(I) + W(I)
20 CONTINUE
CALL SETV(Z, 0.0D0, 0.0D0, 0.0D0)
CALL MPLACE(KSIV, S7AT, PS, Z)
RETURN
END
C
C S7ATT: the stack's attitude, LM (and S-IVB) body axes in S7AT.
C The S-IVB held "a fixed inertial attitude to provide a stable
C docking platform" (MPR-SAT-FE-69-9, printed p. 11-1), reached by
C a manoeuvre that was to be "completed at plus 09 plus 20", so
C "the Sun will shine across the top of the LM after separation"
C (Apollo 11 Flight Journal, 002:54:09 and commentary). The
C attitude itself is our guess: the stack's X axis square to the
C Sun and in the plane of the trajectory, pointing along the
C motion, frozen at 3:09:20; the LM front (+Z) toward the Sun.
SUBROUTINE S7ATT
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T, R(3), V(3), S(3), H(3), X(3), Y(3), VDOT
DOUBLE PRECISION EVGET
INTEGER I
C The attitude time from the scenario (TDATT).
T = EVGET(KETDA)
CALL VSTATE(T, 1, R, V)
CALL SUNG(T, S)
CALL VCRS(R, V, H)
CALL VUNIT(H)
CALL VCRS(S, H, X)
CALL VUNIT(X)
IF (VDOT(X, V) .LT. 0.0D0) CALL SETV(X, -X(1), -X(2), -X(3))
CALL VCRS(S, X, Y)
DO 10 I = 1, 3
S7AT(I,1) = X(I)
S7AT(I,2) = Y(I)
S7AT(I,3) = S(I)
10 CONTINUE
RETURN
END
C
C VSETIN: the view inputs for the next frame (the chassis calls it
C before VFRAME): view, camera target, label level. Out-of-range
C values fall back to 0.
SUBROUTINE VSETIN(IV, IT, IL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER IV, IT, IL
IVIEW = IV
ITARG = IT
ILABL = IL
IF (IVIEW .LT. 0 .OR. IVIEW .GT. 3) IVIEW = 0
IF (ITARG .LT. 0 .OR. ITARG .GT. 5) ITARG = 0
IF (ILABL .LT. 0 .OR. ILABL .GT. 3) ILABL = 0
RETURN
END
C
C-----------------------------------------------------------------------
C S8ATT: the docked stack's attitude in passive thermal control,
C CSM body axes in S8AT. "In this attitude the spacecraft will be
C rotated at a rate of about 3 revolutions per hour" (Apollo 11
C Flight Journal, commentary after 008:11:00); "rotated about its
C X axis" (same, Apollo Control at 8 hours 59 minutes); "PTC is
C started now" (Collins, 010:58:19). The axis is ours: square to
C the ecliptic, so the Sun stays square to it (J2000 ecliptic pole);
C the roll starts with +Z toward the Sun. Before PTC, no roll.
C-----------------------------------------------------------------------
SUBROUTINE S8ATT(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, X(3), Y(3), Z(3), EP, PH, C, S, EVGET
INTEGER I
EP = 23.4392911D0 * DR
CALL SETV(X, 0.0D0, -DSIN(EP), DCOS(EP))
CALL VCRS(X, SUNU, Y)
CALL VUNIT(Y)
CALL VCRS(X, Y, Z)
PH = 0.0D0
IF (GET .GT. EVGET(KEPTC)) PH = 3.0D0 * 2.0D0 * PI / 3600.0D0
& * (GET - EVGET(KEPTC))
C = DCOS(PH)
S = DSIN(PH)
DO 10 I = 1, 3
S8AT(I,1) = X(I)
S8AT(I,2) = C * Y(I) + S * Z(I)
S8AT(I,3) = C * Z(I) - S * Y(I)
10 CONTINUE
RETURN
END
C
C S8POSE: place the docked stack (STKPL), the LM with its gear
C stowed (KLMS) as in translunar coast (press kit p. 103), 60 m
C along the reference boresight from the camera.
SUBROUTINE S8POSE
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION P(3)
INTEGER I
DO 10 I = 1, 3
P(I) = 0.060D0 * BREF(I)
10 CONTINUE
CALL STKPL(S8AT, P, KLMS)
RETURN
END
C
C-----------------------------------------------------------------------
C S9REF: scene 9's reference attitude, the forward horizon view of
C scene 1 turned to the framing of the photograph AS08-14-2383 as
C it is usually shown (the film frame turned a quarter turn
C clockwise, the lunar horizon at the bottom): yaw, pitch, roll
C offsets in LOOK's sense. The three angles are FITTED by us (not
C sourced): they put the Earth's centre where the photograph has
C it and the horizon at its tilt, measured on the ASU scan; the
C Earth's size, its height above the horizon and its phase are
C then checks, not inputs.
C-----------------------------------------------------------------------
SUBROUTINE S9REF
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION YAW, PIT, ROL
INTEGER I
C Fitted 2026-09-30 (see CLAUDE.md, scene 9): the Earth's centre at
C plot (0.8550, -0.4555) deg and the horizon falling 6.634 deg to
C the right, as measured on the scan.
DATA YAW, PIT, ROL / -0.9043D0, 3.7599D0, -6.7873D0 /
CALL VCRS(BREF, UREF, RREF)
CALL VUNIT(RREF)
CALL LOOK(YAW, PIT, ROL)
DO 10 I = 1, 3
BREF(I) = CB(I)
UREF(I) = CU(I)
RREF(I) = CR(I)
10 CONTINUE
RETURN
END
src/ephem.f
C=======================================================================
C
C V I E W - 1 1 0 8 EPHEMERIS
C
C Core element. Time, Sun, Moon and their orientation. One
C relocatable element of the kernel; see vdrive.f for the list.
C
C WORLD MODEL
C Moon: Meeus ch. 47 (ELP-2000/82 abridged, 60 + 60 terms),
C within about 15 km of JPL Horizons over Apollo 8 and 11
C (docs/simulation.md); precessed to J2000 (PRECM). Sun: the
C low-precision series of the Astronomical Almanac, ecliptic of
C date carried to the J2000 equinox by the general precession in
C longitude. Moon orientation: IAU (Archinal et al.) with its
C periodic terms. Earth rotation: GMST about the pole of date,
C carried to J2000 by the IAU 1976 precession (PRECM).
C
C=======================================================================
C
C=======================================================================
C TIME AND EPHEMERIDES
C=======================================================================
SUBROUTINE TSET(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET
C RESTOMOD BEGIN: J2000 epoch (IAU 1976/1984), GMST of Aoki 1982
TJD = TJD0 + GET / 86400.0D0
TDAY = TJD - 2451545.0D0
TCEN = TDAY / 36525.0D0
GMST = DMOD(280.46061837D0 + 360.98564736629D0 * TDAY, 360.0D0)
GMST = GMST * DR
C RESTOMOD END
RETURN
END
C
C-----------------------------------------------------------------------
C MOONG: geocentric Moon, EQ km.
C
C Meeus, Astronomical Algorithms (2nd ed., 1998), ch. 47: the
C ELP-2000/82 lunar theory (Chapront-Touze and Chapront, 1982)
C abridged to the 60 + 60 periodic terms of tables 47.A and 47.B,
C with the E factor for terms in M and the additive terms A1-A3.
C The tables are data (data/meeus47.txt, BLOCK DATA /CMEEUS/).
C Mean ecliptic and equinox of date to J2000 equatorial: mean
C obliquity of date, then the IAU 1976 precession (Lieske 1977;
C Meeus eqs. 22.2 and 21.3, PRECM). Time is TT: UTC plus 32.184 s
C plus
C TAI-UTC by the USNO formula for 1968-02-01 to 1972-01-01,
C 4.21317 s + (MJD - 39126) x 0.002592 s (maia.usno.navy.mil,
C ser7/tai-utc.dat); ours outside that span too.
C Checked against JPL Horizons, docs/simulation.md. MSC's RTCC
C did not compute its ephemeris from a short series: its
C "ephemeris subroutines used in the RTCC will be system
C subroutines", reading "an ephemeris tape" (Analytical Mechanics
C Associates Report 68-4, NAS 9-4036, April 1968, p. 17),
C so a precise Moon is nearer period practice than the
C low-precision formula this replaces.
C-----------------------------------------------------------------------
SUBROUTINE MOONG(GET, P)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, P(3)
DOUBLE PRECISION T, TJ, XMJD, LP, DD, AM, AMP, F, A1, A2, A3, E
DOUBLE PRECISION SL, SR, SB, ARG, EK, L, B, R, EPS, CE, SE
DOUBLE PRECISION X, Y, Z, Q(3), PR(3,3), SND, CSD
INTEGER K, J, M
C RESTOMOD BEGIN: ELP-2000/82 (1982) as abridged by Meeus (1991,
C 1998); IAU 1976 precession; TT and TAI-UTC (1970s definitions)
XMJD = TJD0 + GET / 86400.0D0 - 2400000.5D0
TJ = TJD0 + (GET + 32.184D0 + 4.21317D0
& + (XMJD - 39126.0D0) * 0.002592D0) / 86400.0D0
T = (TJ - 2451545.0D0) / 36525.0D0
LP = 218.3164477D0 + 481267.88123421D0 * T
& - 0.0015786D0 * T**2 + T**3 / 538841.0D0 - T**4 / 65194000.0D0
DD = 297.8501921D0 + 445267.1114034D0 * T
& - 0.0018819D0 * T**2 + T**3 / 545868.0D0
& - T**4 / 113065000.0D0
AM = 357.5291092D0 + 35999.0502909D0 * T
& - 0.0001536D0 * T**2 + T**3 / 24490000.0D0
AMP = 134.9633964D0 + 477198.8675055D0 * T
& + 0.0087414D0 * T**2 + T**3 / 69699.0D0
& - T**4 / 14712000.0D0
F = 93.2720950D0 + 483202.0175233D0 * T
& - 0.0036539D0 * T**2 - T**3 / 3526000.0D0
& + T**4 / 863310000.0D0
A1 = 119.75D0 + 131.849D0 * T
A2 = 53.09D0 + 479264.290D0 * T
A3 = 313.45D0 + 481266.484D0 * T
E = 1.0D0 - 0.002516D0 * T - 0.0000074D0 * T**2
LP = DMOD(LP, 360.0D0)
DD = DMOD(DD, 360.0D0)
AM = DMOD(AM, 360.0D0)
AMP = DMOD(AMP, 360.0D0)
F = DMOD(F, 360.0D0)
SL = 0.0D0
SR = 0.0D0
SB = 0.0D0
DO 10 K = 1, 60
J = 6 * (K - 1)
M = IABS(MMA(J + 2))
EK = 1.0D0
IF (M .EQ. 1) EK = E
IF (M .EQ. 2) EK = E * E
ARG = DBLE(MMA(J + 1)) * DD + DBLE(MMA(J + 2)) * AM
& + DBLE(MMA(J + 3)) * AMP + DBLE(MMA(J + 4)) * F
SL = SL + DBLE(MMA(J + 5)) * EK * SND(ARG)
SR = SR + DBLE(MMA(J + 6)) * EK * CSD(ARG)
10 CONTINUE
DO 20 K = 1, 60
J = 5 * (K - 1)
M = IABS(MMB(J + 2))
EK = 1.0D0
IF (M .EQ. 1) EK = E
IF (M .EQ. 2) EK = E * E
ARG = DBLE(MMB(J + 1)) * DD + DBLE(MMB(J + 2)) * AM
& + DBLE(MMB(J + 3)) * AMP + DBLE(MMB(J + 4)) * F
SB = SB + DBLE(MMB(J + 5)) * EK * SND(ARG)
20 CONTINUE
SL = SL + 3958.0D0 * SND(A1) + 1962.0D0 * SND(LP - F)
& + 318.0D0 * SND(A2)
SB = SB - 2235.0D0 * SND(LP) + 382.0D0 * SND(A3)
& + 175.0D0 * SND(A1 - F) + 175.0D0 * SND(A1 + F)
& + 127.0D0 * SND(LP - AMP) - 115.0D0 * SND(LP + AMP)
L = (LP + SL * 1.0D-6) * DR
B = SB * 1.0D-6 * DR
R = 385000.56D0 + SR * 1.0D-3
C Ecliptic of date to equator of date, mean obliquity (eq. 22.2).
EPS = (23.4392911D0 - 0.0130041667D0 * T
& - 1.6388889D-7 * T**2 + 5.0361111D-7 * T**3) * DR
CE = DCOS(EPS)
SE = DSIN(EPS)
X = R * DCOS(B) * DCOS(L)
Y = R * DCOS(B) * DSIN(L)
Z = R * DSIN(B)
Q(1) = X
Q(2) = CE * Y - SE * Z
Q(3) = SE * Y + CE * Z
C Equator and equinox of date to J2000: PRECM's transpose.
CALL PRECM(T, PR)
CALL MTXV(PR, Q, P)
C RESTOMOD END
RETURN
END
C
C-----------------------------------------------------------------------
C SND: sine of an angle in degrees, reduced first.
C-----------------------------------------------------------------------
DOUBLE PRECISION FUNCTION SND(A)
DOUBLE PRECISION A
DOUBLE PRECISION PI, DR
C RESTOMOD: real-valued PARAMETER; FORTRAN V PARAMETER was
C integer-only (UP-4046 Rev 3, sec. 10.4.1, p. 10-8)
PARAMETER (PI=3.141592653589793D0, DR=PI/180.0D0)
SND = DSIN(DMOD(A, 360.0D0) * DR)
RETURN
END
C
C-----------------------------------------------------------------------
C SUNG: unit vector Earth to Sun, EQ.
C-----------------------------------------------------------------------
SUBROUTINE SUNG(GET, U)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, U(3)
DOUBLE PRECISION D, T, G, L, CE, SE, SND
C RESTOMOD BEGIN: Astronomical Almanac low-precision Sun, 1980s
D = TJD0 + GET / 86400.0D0 - 2451545.0D0
T = D / 36525.0D0
G = 357.528D0 + 0.9856003D0 * D
L = 280.460D0 + 0.9856474D0 * D + 1.915D0 * SND(G)
& + 0.020D0 * SND(2.0D0 * G)
L = DMOD(L - 1.3969713D0 * T, 360.0D0) * DR
CE = DCOS(23.4392911D0 * DR)
SE = DSIN(23.4392911D0 * DR)
U(1) = DCOS(L)
U(2) = CE * DSIN(L)
U(3) = SE * DSIN(L)
C RESTOMOD END
RETURN
END
C
C-----------------------------------------------------------------------
C MOONRT: Moon-fixed to EQ rotation, IAU pole and prime meridian
C with the periodic terms E1..E13. M = RZ(A0+90) RX(90-D0) RZ(W).
C-----------------------------------------------------------------------
SUBROUTINE MOONRT(GET, M)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, M(3,3)
DOUBLE PRECISION D, T, E(13), A0, D0, W, SND, CSD
DOUBLE PRECISION R1(3,3), R2(3,3), R3(3,3), R4(3,3)
C RESTOMOD BEGIN: IAU WGCCRE lunar orientation, 1980s-2010s
D = TJD0 + GET / 86400.0D0 - 2451545.0D0
T = D / 36525.0D0
E(1) = 125.045D0 - 0.0529921D0 * D
E(2) = 250.089D0 - 0.1059842D0 * D
E(3) = 260.008D0 + 13.0120009D0 * D
E(4) = 176.625D0 + 13.3407154D0 * D
E(5) = 357.529D0 + 0.9856003D0 * D
E(6) = 311.589D0 + 26.4057084D0 * D
E(7) = 134.963D0 + 13.0649930D0 * D
E(8) = 276.617D0 + 0.3287146D0 * D
E(9) = 34.226D0 + 1.7484877D0 * D
E(10) = 15.134D0 - 0.1589763D0 * D
E(11) = 119.743D0 + 0.0036096D0 * D
E(12) = 239.961D0 + 0.1643573D0 * D
E(13) = 25.053D0 + 12.9590088D0 * D
A0 = 269.9949D0 + 0.0031D0 * T - 3.8787D0 * SND(E(1))
& - 0.1204D0 * SND(E(2)) + 0.0700D0 * SND(E(3))
& - 0.0172D0 * SND(E(4)) + 0.0072D0 * SND(E(6))
& - 0.0052D0 * SND(E(10)) + 0.0043D0 * SND(E(13))
D0 = 66.5392D0 + 0.0130D0 * T + 1.5419D0 * CSD(E(1))
& + 0.0239D0 * CSD(E(2)) - 0.0278D0 * CSD(E(3))
& + 0.0068D0 * CSD(E(4)) - 0.0029D0 * CSD(E(6))
& + 0.0009D0 * CSD(E(7)) + 0.0008D0 * CSD(E(10))
& - 0.0009D0 * CSD(E(13))
W = 38.3213D0 + 13.17635815D0 * D - 1.4D-12 * D * D
& + 3.5610D0 * SND(E(1)) + 0.1208D0 * SND(E(2))
& - 0.0642D0 * SND(E(3)) + 0.0158D0 * SND(E(4))
& + 0.0252D0 * SND(E(5)) - 0.0066D0 * SND(E(6))
& - 0.0047D0 * SND(E(7)) - 0.0046D0 * SND(E(8))
& + 0.0028D0 * SND(E(9)) + 0.0052D0 * SND(E(10))
& + 0.0040D0 * SND(E(11)) + 0.0019D0 * SND(E(12))
& - 0.0044D0 * SND(E(13))
CALL ROTZ((A0 + 90.0D0) * DR, R1)
CALL ROTX((90.0D0 - D0) * DR, R2)
CALL ROTZ(DMOD(W, 360.0D0) * DR, R3)
CALL MXM(R1, R2, R4)
CALL MXM(R4, R3, M)
C RESTOMOD END
RETURN
END
C
DOUBLE PRECISION FUNCTION CSD(A)
DOUBLE PRECISION A
DOUBLE PRECISION PI, DR
C RESTOMOD: real-valued PARAMETER; FORTRAN V PARAMETER was
C integer-only (UP-4046 Rev 3, sec. 10.4.1, p. 10-8)
PARAMETER (PI=3.141592653589793D0, DR=PI/180.0D0)
CSD = DCOS(DMOD(A, 360.0D0) * DR)
RETURN
END
C
C GMSTAT: Greenwich mean sidereal time (rad) at GET of the current
C scenario; the formula of TSET.
DOUBLE PRECISION FUNCTION GMSTAT(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET
C RESTOMOD BEGIN: GMST of Aoki 1982
GMSTAT = DMOD(280.46061837D0 + 360.98564736629D0 *
& (TJD0 + GET / 86400.0D0 - 2451545.0D0), 360.0D0) * DR
C RESTOMOD END
RETURN
END
C
C MOONV: the Moon's geocentric position P (km) and velocity V
C (km/s) at GET, the velocity by a central difference over 2 s.
SUBROUTINE MOONV(GET, P, V)
DOUBLE PRECISION GET, P(3), V(3), A(3), B(3)
INTEGER I
CALL MOONG(GET, P)
CALL MOONG(GET - 1.0D0, A)
CALL MOONG(GET + 1.0D0, B)
DO 10 I = 1, 3
V(I) = 0.5D0 * (B(I) - A(I))
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C PRECM: the IAU 1976 precession matrix P for T Julian centuries
C from J2000 (Lieske 1977; Meeus, Astronomical Algorithms, 2nd ed.,
C eq. 21.3): P carries J2000 equatorial vectors to the mean
C equator and equinox of date, and its transpose (MTXV) back.
C The Earth-fixed frame turns about the pole of date by GMST, so
C Earth-fixed positions (coastlines, the reports' latitudes and
C longitudes) need the transpose to meet our J2000 stars and Moon:
C about 0.43 deg of precession between 1969 and 2000.
C-----------------------------------------------------------------------
SUBROUTINE PRECM(T, P)
DOUBLE PRECISION T, P(3,3), ZT, ZZ, TH, C1, S1, C2, S2, C3, S3
DOUBLE PRECISION AS
C RESTOMOD BEGIN: IAU 1976 precession (Lieske 1977)
AS = 3.141592653589793D0 / 180.0D0 / 3600.0D0
ZT = (2306.2181D0 * T + 0.30188D0 * T**2 + 0.017998D0 * T**3) * AS
ZZ = (2306.2181D0 * T + 1.09468D0 * T**2 + 0.018203D0 * T**3) * AS
TH = (2004.3109D0 * T - 0.42665D0 * T**2 - 0.041833D0 * T**3) * AS
C RESTOMOD END
C1 = DCOS(ZT)
S1 = DSIN(ZT)
C2 = DCOS(ZZ)
S2 = DSIN(ZZ)
C3 = DCOS(TH)
S3 = DSIN(TH)
P(1,1) = C1 * C3 * C2 - S1 * S2
P(1,2) = -S1 * C3 * C2 - C1 * S2
P(1,3) = -S3 * C2
P(2,1) = C1 * C3 * S2 + S1 * C2
P(2,2) = -S1 * C3 * S2 + C1 * C2
P(2,3) = -S3 * S2
P(3,1) = C1 * S3
P(3,2) = -S1 * S3
P(3,3) = C3
RETURN
END
src/lcoas.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 7 COAS RETICLE
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C S7COAS: the COAS cross hairs, fixed to the CSM. TN D-6853
C (printed p. 12): "Command module and LM windows and optics
C outlines can be simulated." The reticle's pattern and size here
C are our guess: a cross +-6 deg, open in the middle so the target
C shows.
SUBROUTINE S7COAS(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
IVMODE = 0
ISTYLE = 1
CALL OVLINE(VB, NV, -6.0D0, 0.0D0, -1.0D0, 0.0D0)
CALL OVLINE(VB, NV, 1.0D0, 0.0D0, 6.0D0, 0.0D0)
CALL OVLINE(VB, NV, 0.0D0, -6.0D0, 0.0D0, -1.0D0)
CALL OVLINE(VB, NV, 0.0D0, 1.0D0, 0.0D0, 6.0D0)
RETURN
END
src/learth.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 5 EARTH
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C EARTH. Limb, coastlines (Natural Earth, turned by GMST), the
C terminator, and the night side hatched with meridians every
C 10 deg. Everything behind the Moon is dropped.
C=======================================================================
SUBROUTINE DEARTH(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION D, AE, C(3), U(3), P(3), Q(3), X, Y, VNRM
DOUBLE PRECISION OCCL
INTEGER I, K, J, IOK, IP, LMOCC
D = VNRM(EPOS)
IF (D .LE. RE * 1.0001D0) RETURN
AE = DASIN(RE / D)
DO 10 I = 1, 3
U(I) = EPOS(I) / D
C(I) = -U(I)
10 CONTINUE
IF (U(1)*CB(1) + U(2)*CB(2) + U(3)*CB(3) .LT.
& DCOS(DMIN1(THVIEW * DR + AE, PI))) RETURN
C
C Limb.
IVMODE = 2
CALL CIRCLE(VB, NV, EPOS, RE, C, DACOS(RE / D), 360)
C
C Coastlines.
IVMODE = 1
DO 30 K = 1, NCST
IP = 0
DO 20 J = KCST(K), KCST(K+1) - 1
CALL MXV(MEF, CEV(1,J), Q)
DO 15 I = 1, 3
P(I) = EPOS(I) + RE * Q(I)
15 CONTINUE
CALL PEN(VB, NV, P, IP)
IP = 1
20 CONTINUE
30 CONTINUE
C
C Terminator.
CALL CIRCLE(VB, NV, EPOS, RE, SUNU, 0.5D0 * PI, 180)
C
C Night side shading (SHADE), not from low orbit, where the film
C shows none.
IF (ISCN .NE. 3) CALL SHADE(VB, NV, EPOS, RE, 4)
IVMODE = 0
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (MOD(IFLG, 2) .EQ. 1 .AND. AE .LT. 0.3D0 * FOVH * DR) THEN
CALL PROJ(EPOS, X, Y, IOK)
IF (IOK .EQ. 1) THEN
IF (OCCL(EPOS, MPOS, RM) .LE. 0.0D0) THEN
IF (NACT .EQ. 0 .OR. LMOCC(EPOS, 0) .EQ. 0)
& CALL LABEL(LB, NL, X, Y, 4, 0)
END IF
END IF
END IF
C RESTOMOD END
C The launch pad, with a label level set, once the disc is too big
C to carry the EARTH name (DPAD; the rule is ours).
IF (ILABL .GE. 1 .AND. AE .GE. 0.3D0 * FOVH * DR)
& CALL DPAD(VB, NV, LB, NL)
RETURN
END
C
C-----------------------------------------------------------------------
C DPAD: the scenario's launch pad (its PAD card) as a small boxed
C X, the mark the Moon view gives the Apollo 11 landing site
C (DMOON6), and its name (LB kind 9) beside it. A modern
C addition (ours), drawn only with a label level set (in_lablv
C 1-3; DEARTH). The coastlines are geodetic latitudes placed on
C a sphere and turned with the Earth (MEF, VFRAME), so the pad
C goes on the same footing: a geocentric latitude is made
C geodetic first, TAN(GD) = TAN(GC) / (1 - F)**2, F the flattening
C STATEV uses. Hidden on the Earth's far side and behind the
C Moon or a placed model (ISVIS mode 1).
C-----------------------------------------------------------------------
SUBROUTINE DPAD(VB, NV, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), LB(4,MAXL)
INTEGER NV, NL
DOUBLE PRECISION F, FI, LA, U(3), Q(3), P(3), X, Y, W
INTEGER I, IOK, ISVIS
IF (PADCH(8 * (ISN - 1) + 1) .EQ. 0) RETURN
F = 1.0D0 / 298.257D0
FI = SNPLA(ISN) * DR
IF (SNPGC(ISN) .EQ. 1) FI = DATAN(DTAN(FI) / (1.0D0 - F)**2)
LA = SNPLO(ISN) * DR
U(1) = DCOS(FI) * DCOS(LA)
U(2) = DCOS(FI) * DSIN(LA)
U(3) = DSIN(FI)
CALL MXV(MEF, U, Q)
DO 10 I = 1, 3
P(I) = EPOS(I) + RE * Q(I)
10 CONTINUE
IF (P(1)*CB(1) + P(2)*CB(2) + P(3)*CB(3) .LE. 0.0D0) RETURN
IVMODE = 1
IF (ISVIS(P) .EQ. 0) GO TO 90
CALL PROJ(P, X, Y, IOK)
IF (IOK .EQ. 0) GO TO 90
W = 0.012D0 * FOVH
CALL BOXX(VB, NV, X, Y, W)
IF (ILEV .GE. 1) CALL LABEL(LB, NL, X + W, Y + W, 9, 0)
90 IVMODE = 0
RETURN
END
src/lframe.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 1 PLOT FRAME
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C PLOT FRAME. Box at the field edge, ticks inward every 5 deg
C up to a 25 deg field, 10 deg up to 60, else 20 deg, at
C multiples of the step from 0 (the page letters those values).
C=======================================================================
SUBROUTINE DFRAME(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION ST, TL, B, X
INTEGER K, N
C Only when IFLG bit 1 (frame and ticks) is set.
IF (MOD(IFLG / 2, 2) .NE. 1) RETURN
B = BOXH
ST = 20.0D0
IF (2.0D0 * B .LE. 60.0D0) ST = 10.0D0
IF (2.0D0 * B .LE. 25.0D0) ST = 5.0D0
TL = 0.02D0 * B
CALL EMIT(VB, NV, -B, -B, B, -B)
CALL EMIT(VB, NV, B, -B, B, B)
CALL EMIT(VB, NV, B, B, -B, B)
CALL EMIT(VB, NV, -B, B, -B, -B)
N = INT(B / ST + 1.0D-9)
DO 30 K = -N, N
X = DBLE(K) * ST
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (DABS(X) .LT. B - 1.0D-9) THEN
CALL EMIT(VB, NV, X, -B, X, -B + TL)
CALL EMIT(VB, NV, X, B, X, B - TL)
CALL EMIT(VB, NV, -B, X, -B + TL, X)
CALL EMIT(VB, NV, B, X, B - TL, X)
END IF
C RESTOMOD END
30 CONTINUE
RETURN
END
src/llpd.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 9 LPD AND LM WINDOW
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C LM FRONT WINDOW OVERLAY. Landing point designator scale and the
C commander's window frame, fixed to the LM (drawn in reference
C plot degrees, so they move with free-look). TN D-6853 (printed
C p. 7) says the LPD and scribe marks came from LM window
C engineering drawings; we do not have them. The geometry below
C is our reading of the film's descent frames (film seconds 27-35,
C reference/video_frames/descent_t*.png), taking 60 px per 10 deg
C from the edge numbers and Y = 0 at the "0" labels:
C scale line X = 0 from Y = +36 to the sill, small marks every
C 1.125 deg on alternate sides, a longer mark to the
C right every 9 deg down to the lower cross bar;
C cross bars upper at Y = +29.3, +-11 deg, ticks up at 5, 10;
C lower at Y = -18.8, +-6.5 deg, end ticks down;
C window sill near Y = -35, right edge two lines to the
C frame top (the film's frame is cut at Y = +41).
C=======================================================================
SUBROUTINE OVLPD(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION Y, WX(9), WY(9), YU, YL
INTEGER K, J
IVMODE = 0
ISTYLE = 1
YU = 29.3D0
YL = -18.8D0
CALL OVLINE(VB, NV, 0.0D0, 36.0D0, 0.0D0, -35.5D0)
DO 10 J = -14, 48
Y = YL + 1.125D0 * DBLE(J)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (MOD(J + 16, 8) .EQ. 0 .AND. J .GT. 0) THEN
IF (Y .LT. YU) CALL OVLINE(VB, NV, 0.0D0, Y, 2.5D0, Y)
GO TO 10
END IF
IF (MOD(J + 16, 2) .EQ. 0) THEN
CALL OVLINE(VB, NV, -0.6D0, Y, 0.0D0, Y)
ELSE
CALL OVLINE(VB, NV, 0.0D0, Y, 0.6D0, Y)
END IF
C RESTOMOD END
10 CONTINUE
C Upper cross bar, ticks up at 5 and 10 deg each side and the ends.
CALL OVLINE(VB, NV, -11.0D0, YU, 11.0D0, YU)
DO 20 K = -2, 2
IF (K .NE. 0) CALL OVLINE(VB, NV, 5.0D0 * DBLE(K), YU,
& 5.0D0 * DBLE(K), YU + 1.3D0)
20 CONTINUE
C Lower cross bar with end ticks down.
CALL OVLINE(VB, NV, -6.5D0, YL, 6.5D0, YL)
CALL OVLINE(VB, NV, -6.5D0, YL, -6.5D0, YL - 1.5D0)
CALL OVLINE(VB, NV, 6.5D0, YL, 6.5D0, YL - 1.5D0)
C Window frame: sill, then the right-hand edge as two lines.
WX(1) = -60.0D0
WY(1) = -34.8D0
WX(2) = -30.0D0
WY(2) = -34.5D0
WX(3) = -12.0D0
WY(3) = -35.3D0
WX(4) = -1.0D0
WY(4) = -35.7D0
WX(5) = 13.4D0
WY(5) = 38.0D0
WX(6) = 10.0D0
WY(6) = 40.5D0
DO 40 K = 1, 5
CALL OVLINE(VB, NV, WX(K), WY(K), WX(K+1), WY(K+1))
40 CONTINUE
WX(1) = -1.0D0
WY(1) = -35.7D0
WX(2) = 2.0D0
WY(2) = -33.5D0
WX(3) = 5.5D0
WY(3) = -27.0D0
WX(4) = 16.0D0
WY(4) = 37.0D0
WX(5) = 15.5D0
WY(5) = 38.5D0
WX(6) = 10.0D0
WY(6) = 40.5D0
DO 50 K = 1, 5
CALL OVLINE(VB, NV, WX(K), WY(K), WX(K+1), WY(K+1))
50 CONTINUE
RETURN
END
src/lmoon.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 4 MOON AND CRATERS
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C MOON. Limb (the horizon when close), gazetteer craters, and
C near the surface seeded small craters.
C=======================================================================
SUBROUTINE DMOON(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION D, AM, C(3), U(3), H, HOR, X, Y, OCCL, VNRM
DOUBLE PRECISION G(3), Q(3), CA
INTEGER I, IOK, K
D = VNRM(MPOS)
IF (D .LE. RM * 1.000001D0) RETURN
AM = DASIN(RM / D)
DO 10 I = 1, 3
U(I) = MPOS(I) / D
C(I) = -U(I)
10 CONTINUE
IF (U(1)*CB(1) + U(2)*CB(2) + U(3)*CB(3) .LT.
& DCOS(DMIN1(THVIEW * DR + AM, PI))) RETURN
H = D - RM
HOR = DACOS(RM / D)
C
C Limb.
IVMODE = 5
CALL CIRCLE(VB, NV, MPOS, RM, C, HOR, 720)
C
C Gazetteer craters inside the visible cap.
IVMODE = 3
CA = DCOS(DMIN1(HOR + 0.05D0, PI))
DO 20 K = 1, NCRAT
IF (CRV(1,K)*CAMF(1) + CRV(2,K)*CAMF(2) + CRV(3,K)*CAMF(3)
& .LT. CA * D) GO TO 20
C Whole-disc view: gazetteer craters of 25 km and up only (our
C floor, so the disc is not a solid mass), unlabelled.
IF (ISCN .EQ. 6 .AND. CRDIA(K) .LT. 25.0D0) GO TO 20
CALL CRATER(VB, NV, CRV(1,K), 0.5D0 * CRDIA(K) / RM, IOK)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IOK .EQ. 1 .AND. CRDIA(K) .GE. 20.0D0 .AND.
& MOD(IFLG, 2) .EQ. 1 .AND. ISCN .NE. 6) THEN
DO 15 I = 1, 3
G(I) = RM * CRV(I,K)
15 CONTINUE
CALL MXV(MMF, G, Q)
DO 16 I = 1, 3
Q(I) = Q(I) + MPOS(I)
16 CONTINUE
CALL PROJ(Q, X, Y, IOK)
CALL LABEL(LB, NL, X, Y, 2, K)
END IF
C RESTOMOD END
20 CONTINUE
C
C Seeded craters, a FIXED set on the ground (no level of detail
C by range). Sources: VIEW had "two-dimensional crater models"
C built from photographs of the area near Apollo landing site 2
C (TN D-6853, printed p. 7), and "the smallest craters depicted
C ... have a size of 1 minute of arc (1658 ft)" (MSC IN
C 69-FM-197, sec. 3.1). Our model of that, not the original:
C regional patch, site 2 (see DMOON6) +-5 deg, cells of
C 0.1 deg, craters 1658 ft (0.505 km) to 4 km, where the
C gazetteer takes over;
C global fill elsewhere, cells of 0.5 deg, 2.5 to 25 km, so
C orbital views keep the density the film shows (t04).
C The IN (sec. 3.6, 3.7) notes few craters near the site, from
C the flat approach, the oblique look and few large craters;
C the patch density is set low to match. Cells are seeded from
C their indices, so craters stay put from frame to frame.
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (ISCN .EQ. 6) THEN
CALL DMOON6(VB, NV, LB, NL)
ELSE
CALL PCRAT(VB, NV, 0.5D0, 1.2D0, 2.5D0, 25.0D0, 1,
& -90.0D0, 90.0D0, -180.0D0, 180.0D0, 1)
CALL PCRAT(VB, NV, 0.1D0, 0.8D0, 0.505D0, 4.0D0, 2,
& 0.67416D0 - 5.0D0, 0.67416D0 + 5.0D0,
& 23.47314D0 - 5.0D0, 23.47314D0 + 5.0D0, 0)
END IF
C RESTOMOD END
IVMODE = 0
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (MOD(IFLG, 2) .EQ. 1 .AND. AM .LT. 0.3D0 * FOVH * DR) THEN
CALL PROJ(MPOS, X, Y, IOK)
IF (IOK .EQ. 1) THEN
IF (OCCL(MPOS, EPOS, RE) .LE. 0.0D0)
& CALL LABEL(LB, NL, X, Y, 5, 0)
END IF
END IF
C RESTOMOD END
RETURN
END
C
C LLUNIT: selenographic latitude, east longitude (deg) to a unit
C vector (MF). SURFPT: unit MF direction to the surface point,
C camera relative EQ.
SUBROUTINE LLUNIT(FI, LA, U)
DOUBLE PRECISION FI, LA, U(3), DR
DR = 3.141592653589793D0 / 180.0D0
U(1) = DCOS(FI * DR) * DCOS(LA * DR)
U(2) = DCOS(FI * DR) * DSIN(LA * DR)
U(3) = DSIN(FI * DR)
RETURN
END
C
SUBROUTINE SURFPT(CM, P)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION CM(3), P(3), G(3)
INTEGER I
DO 10 I = 1, 3
G(I) = RM * CM(I)
10 CONTINUE
CALL MXV(MMF, G, P)
DO 20 I = 1, 3
P(I) = P(I) + MPOS(I)
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C CRATER: rim circle of angular radius A (rad) about unit centre
C CM (MF). Culled when off frame, beyond the horizon or too small
C to see. IOK = 1 if drawn.
C-----------------------------------------------------------------------
SUBROUTINE CRATER(VB, NV, CM, A, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), CM(3), A
INTEGER NV, IOK
DOUBLE PRECISION D(3), DD, RA, CS, E1(3), E2(3), G(3), P(3)
DOUBLE PRECISION CA, SA, T, CT, ST, SIZ
INTEGER I, K, N
IOK = 0
DO 10 I = 1, 3
D(I) = RM * CM(I) - CAMF(I)
10 CONTINUE
DD = DSQRT(D(1)*D(1) + D(2)*D(2) + D(3)*D(3))
RA = RM * A / DD
C Apparent diameter against the field. Rims under 1.2 percent
C of the frame (about 12 points of a 1024-point recorder raster,
C docs/univac-1108.md) are not drawn: our choice of floor.
SIZ = 2.0D0 * RA / DR / (2.0D0 * FOVH)
IF (SIZ .LT. 0.012D0) RETURN
C Outside the cone through the frame corners (COS(T+RA) is at
C least COS(T) - RA).
CS = (D(1)*CBMF(1) + D(2)*CBMF(2) + D(3)*CBMF(3)) / DD
IF (CS .LT. CSVIEW - RA) RETURN
N = 10 + INT(SIZ * 60.0D0)
IF (N .GT. 36) N = 36
CALL PERP(CM, E1, E2)
CA = DCOS(A) * RM
SA = DSIN(A) * RM
DO 30 K = 0, N
T = DBLE(K) * 2.0D0 * PI / DBLE(N)
CT = DCOS(T) * SA
ST = DSIN(T) * SA
DO 20 I = 1, 3
G(I) = CA * CM(I) + CT * E1(I) + ST * E2(I)
20 CONTINUE
CALL MXV(MMF, G, P)
DO 25 I = 1, 3
P(I) = P(I) + MPOS(I)
25 CONTINUE
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (K .EQ. 0) THEN
CALL PEN(VB, NV, P, 0)
ELSE
CALL PEN(VB, NV, P, 1)
END IF
C RESTOMOD END
30 CONTINUE
IOK = 1
RETURN
END
C
C-----------------------------------------------------------------------
C PCRAT: seeded craters on a latitude-longitude grid of GC deg,
C inside the fixed box F1..F2 lat, L1..L2 east lon (deg); IEXC=1
C skips cells inside the site 2 patch (drawn by its own level).
C AVG mean craters per cell, diameters DMIN..DMAX km with a -2
C power law. Only cells under a 7 x 7 grid of sight lines across
C the frame are visited, and only craters on the near side of the
C horizon are drawn: culling by visibility, not by range. More
C than 40000 cells in view is taken as "Moon too small to show
C seeded craters" (every one would fall under the size floor in
C CRATER). Random numbers: Park-Miller minimal standard
C generator, exact in double precision.
C-----------------------------------------------------------------------
SUBROUTINE PCRAT(VB, NV, GC, AVG, DMIN, DMAX, LEV,
& F1, F2, L1, L2, IEXC)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), GC, AVG, DMIN, DMAX
DOUBLE PRECISION F1, F2, L1, L2, P1, CH, DCAM, FC, LC, CFR
INTEGER NV, LEV, IEXC
DOUBLE PRECISION U(3), W(3), GP(3), B, C, DS, T, FI, DL, LO0
DOUBLE PRECISION FMIN, FMAX, LMIN, LMAX, MARG, SEED, RND
DOUBLE PRECISION CLAT0, CLON0, CF, DIA, CM(3), Q
INTEGER I, J, K, IA, IB, JA, JB, NLON, JJ, NC, IOK, IFLR
C
C The size floor of CRATER applied to the largest crater of this
C level at the nearest ground: if even that is too small, stop.
DCAM = DSQRT(CAMF(1)**2 + CAMF(2)**2 + CAMF(3)**2)
IF (DMAX / (DCAM - RM) / DR / (2.0D0 * FOVH) .LT. 0.012D0) RETURN
LO0 = DATAN2(CAMF(2), CAMF(1)) / DR
FMIN = 90.0D0
FMAX = -90.0D0
LMIN = 180.0D0
LMAX = -180.0D0
C = CAMF(1)**2 + CAMF(2)**2 + CAMF(3)**2 - RM * RM
DO 20 I = 0, 6
DO 10 J = 0, 6
CALL UNPROJ(BOXH * DBLE(I - 3) / 3.0D0,
& BOXH * DBLE(J - 3) / 3.0D0, 0, W)
CALL MTXV(MMF, W, U)
B = U(1) * CAMF(1) + U(2) * CAMF(2) + U(3) * CAMF(3)
DS = B * B - C
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (DS .GE. 0.0D0 .AND. -B - DSQRT(DS) .GT. 0.0D0) THEN
T = -B - DSQRT(DS)
ELSE
T = -B
IF (T .LT. 0.0D0) T = 0.0D0
END IF
C RESTOMOD END
DO 5 K = 1, 3
GP(K) = CAMF(K) + T * U(K)
5 CONTINUE
CALL VUNIT(GP)
FI = DASIN(GP(3)) / DR
DL = DATAN2(GP(2), GP(1)) / DR - LO0
IF (DL .GT. 180.0D0) DL = DL - 360.0D0
IF (DL .LT. -180.0D0) DL = DL + 360.0D0
FMIN = DMIN1(FMIN, FI)
FMAX = DMAX1(FMAX, FI)
LMIN = DMIN1(LMIN, DL)
LMAX = DMAX1(LMAX, DL)
10 CONTINUE
20 CONTINUE
C Take in the ground under the camera too, then pad.
FI = DASIN(CAMF(3) / DSQRT(C + RM * RM)) / DR
FMIN = DMIN1(FMIN, FI)
FMAX = DMAX1(FMAX, FI)
LMIN = DMIN1(LMIN, 0.0D0)
LMAX = DMAX1(LMAX, 0.0D0)
MARG = 0.5D0 * DMAX / RM / DR + GC
FMIN = DMAX1(FMIN - MARG, -89.0D0)
FMAX = DMIN1(FMAX + MARG, 89.0D0)
LMIN = LMIN - MARG / DMAX1(DCOS(DMAX1(DABS(FMIN),
& DABS(FMAX)) * DR), 0.05D0)
LMAX = LMAX + MARG /
& DMAX1(DCOS(DMAX1(DABS(FMIN), DABS(FMAX)) * DR), 0.05D0)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (LMAX - LMIN .GT. 360.0D0) THEN
LMIN = -180.0D0
LMAX = 180.0D0
END IF
C RESTOMOD END
C Clip to the fixed box of this level.
FMIN = DMAX1(FMIN, F1)
FMAX = DMIN1(FMAX, F2)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (L2 - L1 .LT. 360.0D0) THEN
P1 = DMOD(L1 - LO0 + 540.0D0, 360.0D0) - 180.0D0
LMIN = DMAX1(LMIN, P1)
LMAX = DMIN1(LMAX, P1 + (L2 - L1))
END IF
C RESTOMOD END
IF (FMIN .GT. FMAX .OR. LMIN .GT. LMAX) RETURN
DCAM = DSQRT(C + RM * RM)
CH = RM / DCAM
IA = IFLR((FMIN + 90.0D0) / GC)
IB = IFLR((FMAX + 90.0D0) / GC)
JA = IFLR((LO0 + LMIN + 180.0D0) / GC)
JB = IFLR((LO0 + LMAX + 180.0D0) / GC)
IF (DBLE(IB - IA + 1) * DBLE(JB - JA + 1) .GT. 40000.0D0) RETURN
NLON = NINT(360.0D0 / GC)
DO 60 I = IA, IB
C Thin the cells toward the poles, where they shrink in area.
CFR = DCOS((DBLE(I) + 0.5D0) * GC * DR - 0.5D0 * PI)
DO 50 J = JA, JB
JJ = MOD(J, NLON)
IF (JJ .LT. 0) JJ = JJ + NLON
C Seed from both cell indices through a nonlinear mix (the
C fraction of a square), so neighbouring cells do not start the
C linear generator on a lattice (that showed as crater rows).
Q = DBLE(I) * 0.6180339887D0 + DBLE(JJ) * 0.7548776662D0
& + DBLE(LEV) * 0.5698402910D0
Q = DMOD(Q * Q * 7919.0D0, 1.0D0)
C RESTOMOD: seed for the Park-Miller generator (1988)
SEED = 1.0D0 + DBLE(INT(Q * 2147483645.0D0))
Q = RND(SEED)
NC = INT(RND(SEED) * (2.0D0 * AVG + 1.0D0))
CLAT0 = DBLE(I) * GC - 90.0D0
CLON0 = DBLE(JJ) * GC - 180.0D0
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IEXC .EQ. 1) THEN
FC = CLAT0 + 0.5D0 * GC - 0.67416D0
LC = CLON0 + 0.5D0 * GC - 23.47314D0
IF (DABS(FC) .LT. 5.0D0 .AND. DABS(LC) .LT. 5.0D0) GO TO 50
END IF
C RESTOMOD END
DO 40 K = 1, NC
FI = (CLAT0 + RND(SEED) * GC) * DR
DL = (CLON0 + RND(SEED) * GC) * DR
Q = RND(SEED)
CF = RND(SEED)
IF (Q .GT. CFR) GO TO 40
DIA = DMIN / DSQRT(1.0D0 - CF * (1.0D0 - (DMIN/DMAX)**2))
CM(1) = DCOS(FI) * DCOS(DL)
CM(2) = DCOS(FI) * DSIN(DL)
CM(3) = DSIN(FI)
C Near side of the horizon only (centre within the cap
C seen from the camera, padded by the crater's radius).
IF ((CM(1)*CAMF(1) + CM(2)*CAMF(2) + CM(3)*CAMF(3)) / DCAM
& .LT. CH - 0.5D0 * DIA / RM) GO TO 40
CALL CRATER(VB, NV, CM, 0.5D0 * DIA / RM, IOK)
40 CONTINUE
50 CONTINUE
60 CONTINUE
RETURN
END
C
DOUBLE PRECISION FUNCTION RND(SEED)
DOUBLE PRECISION SEED
C RESTOMOD BEGIN: Park-Miller minimal standard generator, 1988
SEED = DMOD(16807.0D0 * SEED, 2147483647.0D0)
RND = SEED / 2147483647.0D0
C RESTOMOD END
RETURN
END
C
INTEGER FUNCTION IFLR(X)
DOUBLE PRECISION X
IFLR = INT(X)
IF (DBLE(IFLR) .GT. X) IFLR = IFLR - 1
RETURN
END
src/lmoon6.f
C=======================================================================
C
C V I E W - 1 1 0 8 MOON VIEW
C
C Part of layer 4: the whole-disc Moon view's extras,
C called by DMOON in scene 6. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C-----------------------------------------------------------------------
C DMOON6: the whole-disc Moon view's extras. Terminator; night
C side shading as for the Earth (TN D-6853, printed p. 8); maria,
C lacus, sinus and oceanus from the IAU gazetteer, which gives
C only a centre and a diameter, so each is drawn as a circle of
C that diameter, not its true outline; the Apollo 11 landing site
C as a small boxed X. Labels (LB kind 6 = mare, id its index;
C kind 7 = landing site) when IFLG bit 0 is set.
C-----------------------------------------------------------------------
SUBROUTINE DMOON6(VB, NV, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), LB(4,MAXL)
INTEGER NV, NL
DOUBLE PRECISION CM(3), P(3), X, Y, W, DCAM, HN, WD
DOUBLE PRECISION XA(64), XB(64), YA(64), YB(64)
INTEGER I, K, IOK, ISVIS, NP, NC
IVMODE = 3
CALL CIRCLE(VB, NV, MPOS, RM, SUNU, 0.5D0 * PI, 360)
CALL SHADE(VB, NV, MPOS, RM, 6)
DCAM = DSQRT(CAMF(1)**2 + CAMF(2)**2 + CAMF(3)**2)
C Labels placed so far, as text rectangles in plot deg (0.7 name
C height per character wide, one name height tall, the extents
C TXALL will give them), for the clutter test: a label whose
C rectangle meets one already placed is dropped. Our rule.
NP = 0
HN = 0.028D0 * FOVH
C Apollo 11 landing site, first so its label always wins. The
C 1969 reports call it "landing site 2" (MSC IN 69-FM-197). LM
C position 0.67416 N, 23.47314 E, planetocentric Mean Earth/Polar
C Axis (DE421), from LRO images: NSSDC, "Apollo Landing Site
C Coordinates", https://nssdc.gsfc.nasa.gov/planetary/lunar/
C lunar_sites.html, citing Wagner et al., Icarus 283, 92-103
C (2017). The same values are in the run deck's SITE card (for
C the descent) and in the crater patch (lmoon.f).
IVMODE = 3
CALL LLUNIT(0.67416D0, 23.47314D0, CM)
CALL SURFPT(CM, P)
IF (ISVIS(P) .EQ. 0) GO TO 10
CALL PROJ(P, X, Y, IOK)
W = 0.012D0 * FOVH
CALL BOXX(VB, NV, X, Y, W)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (MOD(IFLG, 2) .EQ. 1) THEN
CALL LABEL(LB, NL, X + W, Y + W, 7, 0)
NP = 1
XA(1) = X + W + 0.4D0 * HN
XB(1) = XA(1) + 0.7D0 * HN * 22.0D0
YA(1) = Y + W + 0.4D0 * HN
YB(1) = YA(1) + HN
END IF
C RESTOMOD END
C Maria etc., largest first (the table is sorted by diameter).
C All are drawn; only a mare or oceanus of any size, or another
C feature of 150 km or more, is labelled (our clutter rule).
10 DO 30 K = 1, NMARE
CALL LLUNIT(MRLAT(K), MRLON(K), CM)
CALL CRATER(VB, NV, CM, 0.5D0 * MRDIA(K) / RM, IOK)
IF (MOD(IFLG, 2) .EQ. 0) GO TO 30
IF (MRDIA(K) .LT. 150.0D0 .AND. MRCH((K - 1) * 24 + 1) .NE. 77
& .AND. MRCH((K - 1) * 24 + 1) .NE. 79) GO TO 30
IF (CM(1)*CAMF(1) + CM(2)*CAMF(2) + CM(3)*CAMF(3)
& .LT. RM * RM / DCAM) GO TO 30
CALL SURFPT(CM, P)
CALL PROJ(P, X, Y, IOK)
C Name length, then its rectangle, centred on the point.
NC = 0
DO 15 I = 1, 24
IF (MRCH((K - 1) * 24 + I) .EQ. 0) GO TO 16
NC = NC + 1
15 CONTINUE
16 WD = 0.35D0 * HN * DBLE(NC)
DO 20 I = 1, NP
IF (X - WD .LT. XB(I) .AND. X + WD .GT. XA(I) .AND.
& Y - 0.5D0 * HN .LT. YB(I) .AND. Y + 0.5D0 * HN .GT. YA(I))
& GO TO 30
20 CONTINUE
CALL LABEL(LB, NL, X, Y, 6, K)
IF (NP .GE. 64) GO TO 30
NP = NP + 1
XA(NP) = X - WD
XB(NP) = X + WD
YA(NP) = Y - 0.5D0 * HN
YB(NP) = Y + 0.5D0 * HN
30 CONTINUE
IVMODE = 0
RETURN
END
src/lshad.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 8 LM SHADOW
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C-----------------------------------------------------------------------
C LMSHAD: the LM's shadow on the ground in the descent. Every
C vertex of the LM wireframe (model KLMD, MLIB) is carried along
C the Sun's
C direction to the lunar sphere and the edges are drawn there as a
C surface feature (facing test, so it hides below the horizon and
C foreshortens like a crater). The film shows a small LM-shaped
C figure below the horizon late in the descent (descent_t35.png);
C that it is the shadow is our reading. The model's pads (X =
C -1.1 m, LMGEAR) sit at the footpads, LMEYE below the eye.
C-----------------------------------------------------------------------
SUBROUTINE LMSHAD(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION PMF(3), XB(3), YB(3), ZB(3), SMF(3), V(3)
DOUBLE PRECISION A(3), B(3)
INTEGER IS, J, K, IOK
CALL LMDESC(GET, PMF, XB, YB, ZB)
CALL MTXV(MMF, SUNU, SMF)
IVMODE = 3
ISTYLE = 1
DO 20 IS = MDS1(KLMD), MDS2(KLMD)
DO 10 J = 1, NLE(IS)
DO 5 K = 1, 3
V(K) = LMV(K,LME(1,J,IS),IS)
5 CONTINUE
CALL SHADPT(V, PMF, XB, YB, ZB, SMF, A, IOK)
IF (IOK .EQ. 0) GO TO 10
DO 6 K = 1, 3
V(K) = LMV(K,LME(2,J,IS),IS)
6 CONTINUE
CALL SHADPT(V, PMF, XB, YB, ZB, SMF, B, IOK)
IF (IOK .EQ. 0) GO TO 10
CALL PEN(VB, NV, A, 0)
CALL PEN(VB, NV, B, 1)
10 CONTINUE
20 CONTINUE
DO 30 J = MDX1(KLMD), MDX2(KLMD)
IF (LXS(J) .NE. 0) GO TO 30
CALL SHADPT(LXL(1,J), PMF, XB, YB, ZB, SMF, A, IOK)
IF (IOK .EQ. 0) GO TO 30
CALL SHADPT(LXL(4,J), PMF, XB, YB, ZB, SMF, B, IOK)
IF (IOK .EQ. 0) GO TO 30
CALL PEN(VB, NV, A, 0)
CALL PEN(VB, NV, B, 1)
30 CONTINUE
IVMODE = 0
RETURN
END
C
C SHADPT: LM body point V (m; X up, Y right, Z forward) to its
C shadow on the Moon, returned camera relative EQ in P.
SUBROUTINE SHADPT(V, PMF, XB, YB, ZB, SMF, P, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION V(3), PMF(3), XB(3), YB(3), ZB(3), SMF(3), P(3)
DOUBLE PRECISION W(3), G(3), B, C, DS, T
INTEGER IOK, K
IOK = 0
DO 10 K = 1, 3
W(K) = PMF(K) + ((V(1) + 1.1D0 - LMEYE * 1.0D3) * XB(K)
& + V(2) * YB(K)
& + V(3) * ZB(K)) * 1.0D-3
10 CONTINUE
B = W(1) * SMF(1) + W(2) * SMF(2) + W(3) * SMF(3)
C = W(1)**2 + W(2)**2 + W(3)**2 - RM * RM
DS = B * B - C
IF (DS .LT. 0.0D0) RETURN
T = B - DSQRT(DS)
IF (T .LT. 0.0D0) RETURN
DO 20 K = 1, 3
G(K) = W(K) - T * SMF(K)
20 CONTINUE
CALL MXV(MMF, G, P)
DO 30 K = 1, 3
P(K) = P(K) + MPOS(K)
30 CONTINUE
IOK = 1
RETURN
END
src/lstars.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 2 STARS
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C STARS. Points inside the frame not behind the Earth or Moon.
C Nav stars (1..37) labelled when IFLG bit 0 is set.
C=======================================================================
SUBROUTINE DSTARS(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION U(3), X, Y, RAYHIT
INTEGER I, IOK, LMOCC
DO 10 I = 1, NSTAR
U(1) = STX(I)
U(2) = STY(I)
U(3) = STZ(I)
IF (U(1)*CB(1) + U(2)*CB(2) + U(3)*CB(3) .LT. CSVIEW) GO TO 10
CALL PROJ(U, X, Y, IOK)
IF (DABS(X) .GT. BOXH .OR. DABS(Y) .GT. BOXH) GO TO 10
IF (RAYHIT(U, EPOS, RE) .GT. 0.0D0) GO TO 10
IF (RAYHIT(U, MPOS, RM) .GT. 0.0D0) GO TO 10
C Behind a placed spacecraft model (U taken as a point 1 km out).
IF (NACT .EQ. 0) GO TO 8
IF (LMOCC(U, 0) .EQ. 1) GO TO 10
8 IF (NS .GE. MAXS) RETURN
NS = NS + 1
SB(1,NS) = X
SB(2,NS) = Y
SB(3,NS) = STM(I)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (I .LE. NNAV .AND. MOD(IFLG, 2) .EQ. 1) THEN
CALL LABEL(LB, NL, X, Y, 1, I)
END IF
C RESTOMOD END
10 CONTINUE
RETURN
END
src/lsun.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 3 SUN
C
C Layer element. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C SUN. A circle of the Sun's apparent size where it is in view.
C=======================================================================
SUBROUTINE DSUN(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
DOUBLE PRECISION X, Y, RAYHIT, A, C, S
INTEGER IOK, K, LMOCC
IF (SUNU(1)*CB(1) + SUNU(2)*CB(2) + SUNU(3)*CB(3) .LT. CSVIEW)
& RETURN
IF (RAYHIT(SUNU, EPOS, RE) .GT. 0.0D0) RETURN
IF (RAYHIT(SUNU, MPOS, RM) .GT. 0.0D0) RETURN
IF (NACT .EQ. 0) GO TO 5
IF (LMOCC(SUNU, 0) .EQ. 1) RETURN
5 CALL PROJ(SUNU, X, Y, IOK)
IF (IOK .EQ. 0) RETURN
DO 10 K = 0, 23
A = DBLE(K) * PI / 12.0D0
C = DCOS(A) * 0.267D0
S = DSIN(A) * 0.267D0
CALL EMIT(VB, NV, X + C, Y + S,
& X + DCOS(A + PI / 12.0D0) * 0.267D0,
& Y + DSIN(A + PI / 12.0D0) * 0.267D0)
10 CONTINUE
IF (MOD(IFLG, 2) .EQ. 1) CALL LABEL(LB, NL, X, Y, 3, 0)
RETURN
END
src/lvehic.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER 6 VEHICLES
C
C Layer element. Placing and drawing the
C spacecraft models, with hidden lines. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C MCLEAR: no model placed.
SUBROUTINE MCLEAR
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K
DO 10 K = 1, MMOD
MDON(K) = 0
10 CONTINUE
DO 20 K = 1, MSOL
LACT(K) = 0
LINS(K) = 0
20 CONTINUE
NACT = 0
RETURN
END
C
C-----------------------------------------------------------------------
C MPLACE: place model K for this frame. AT: its body axes in EQ
C (columns X, Y, Z). Body point BO (m) goes to P (km, camera
C relative). Its solids go to camera-relative EQ km (LWV, LWN,
C LWD) and join the solids that hide things.
C A model can ride on the observer's own vehicle, a cabin seen
C from inside: AT the vehicle's body axes, BO the eye point in
C body metres, P = 0. A solid with the camera inside it (every
C face turned away) does not hide anything and its edges are all
C drawn (LINS); a cabin is better built of free lines anyway.
C-----------------------------------------------------------------------
SUBROUTINE MPLACE(K, AT, P, BO)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K
DOUBLE PRECISION AT(3,3), P(3), BO(3), V(3), W(3)
INTEGER I, J, IS
DO 10 I = 1, 3
MDP(I,K) = P(I)
MDBO(I,K) = BO(I)
DO 5 J = 1, 3
MDAT(I,J,K) = AT(I,J)
5 CONTINUE
10 CONTINUE
MDON(K) = 1
IF (MDS2(K) .LT. MDS1(K)) RETURN
DO 60 IS = MDS1(K), MDS2(K)
DO 40 J = 1, NLV(IS)
DO 30 I = 1, 3
V(I) = (LMV(I,J,IS) - BO(I)) * 1.0D-3
30 CONTINUE
CALL MXV(AT, V, W)
DO 35 I = 1, 3
LWV(I,J,IS) = P(I) + W(I)
35 CONTINUE
40 CONTINUE
LINS(IS) = 1
DO 50 J = 1, NLF(IS)
CALL MXV(AT, LMN(1,J,IS), W)
DO 45 I = 1, 3
LWN(I,J,IS) = W(I)
45 CONTINUE
LWD(J,IS) = (LMD(J,IS) - LMN(1,J,IS) * BO(1)
& - LMN(2,J,IS) * BO(2) - LMN(3,J,IS) * BO(3)) * 1.0D-3
& + W(1) * P(1) + W(2) * P(2) + W(3) * P(3)
IF (LWD(J,IS) .LE. 0.0D0) LINS(IS) = 0
50 CONTINUE
C An outline model's solids hide nothing (MDHL = 0).
IF (MDHL(K) .EQ. 0) GO TO 60
LACT(IS) = 1
NACT = NACT + 1
60 CONTINUE
RETURN
END
C
C MDRALL: draw every placed model, then their labels and the
C markers of vehicles too small to see (VLABEL).
SUBROUTINE MDRALL(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
INTEGER K
DO 10 K = 1, NMOD
IF (MDON(K) .EQ. 1) CALL MDRAW(VB, NV, K)
10 CONTINUE
ISTYLE = 1
C Vehicle labels and markers, with a label level set (lvlab.f).
IF (ILABL .GE. 1) CALL VLABEL(GET, VB, NV, LB, NL)
RETURN
END
C
C MDRAW: draw placed model K: solid edges, then free lines and
C face marks, each against every placed solid.
SUBROUTINE MDRAW(VB, NV, K)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV)
INTEGER NV, K
DOUBLE PRECISION V(3), W(3), A(3), B(3)
INTEGER I, J, L, IS, IH
C Solid edges. Hidden when both faces are turned away, unless
C the camera is inside the solid.
IF (MDS2(K) .LT. MDS1(K)) GO TO 90
DO 80 IS = MDS1(K), MDS2(K)
DO 70 J = 1, NLE(IS)
DO 65 I = 1, 3
A(I) = LWV(I,LME(1,J,IS),IS)
B(I) = LWV(I,LME(2,J,IS),IS)
65 CONTINUE
IH = 0
IF (LWD(LME(3,J,IS),IS) .GE. 0.0D0 .AND.
& LWD(LME(4,J,IS),IS) .GE. 0.0D0) IH = 1
IF (LINS(IS) .EQ. 1) IH = 0
C An outline model (MDHL = 0): every edge, nothing hidden.
IF (MDHL(K) .EQ. 0) GO TO 68
CALL LMSEG(VB, NV, A, B, IS, IH)
GO TO 70
68 ISTYLE = 1
CALL MSEG(VB, NV, A, B)
70 CONTINUE
80 CONTINUE
C Free lines and face marks.
90 IF (MDX2(K) .LT. MDX1(K)) RETURN
DO 100 J = MDX1(K), MDX2(K)
DO 85 L = 0, 1
DO 82 I = 1, 3
V(I) = (LXL(I + 3 * L, J) - MDBO(I,K)) * 1.0D-3
82 CONTINUE
CALL MXV(MDAT(1,1,K), V, W)
DO 84 I = 1, 3
IF (L .EQ. 0) A(I) = MDP(I,K) + W(I)
IF (L .EQ. 1) B(I) = MDP(I,K) + W(I)
84 CONTINUE
85 CONTINUE
IS = LXS(J)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IS .GT. 0) THEN
IF (LWD(LXF(J),IS) .GE. 0.0D0 .AND. LINS(IS) .EQ. 0)
& GO TO 100
END IF
C RESTOMOD END
IF (MDHL(K) .EQ. 0) GO TO 95
CALL LMSEG(VB, NV, A, B, IS, 0)
GO TO 100
95 ISTYLE = 1
CALL MSEG(VB, NV, A, B)
100 CONTINUE
ISTYLE = 1
RETURN
END
C
C LMSEG: edge A-B of solid IS (0 for a free line) in 12 pieces,
C each tested against the other solids. IHID = 1: the whole edge
C is known hidden.
SUBROUTINE LMSEG(VB, NV, A, B, IS, IHID)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), A(3), B(3)
INTEGER NV, IS, IHID
DOUBLE PRECISION P0(3), P1(3), M(3), F0, F1
INTEGER I, K, NP, IV, IV0, LMOCC
INTEGER IDSH
NP = 12
IDSH = MOD(IFLG / 4, 2)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IHID .EQ. 1) THEN
IF (IDSH .EQ. 1) THEN
ISTYLE = 2
CALL MSEG(VB, NV, A, B)
ISTYLE = 1
END IF
RETURN
END IF
C RESTOMOD END
C Runs of equal visibility are merged into one vector.
IV0 = -1
DO 30 K = 1, NP
F0 = DBLE(K - 1) / DBLE(NP)
F1 = DBLE(K) / DBLE(NP)
DO 10 I = 1, 3
M(I) = A(I) + 0.5D0 * (F0 + F1) * (B(I) - A(I))
10 CONTINUE
IV = 1 - LMOCC(M, IS)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IV .NE. IV0) THEN
IF (IV0 .GE. 0) CALL LMRUN(VB, NV, P0, P1, IV0, IDSH)
DO 15 I = 1, 3
P0(I) = A(I) + F0 * (B(I) - A(I))
15 CONTINUE
IV0 = IV
END IF
C RESTOMOD END
DO 20 I = 1, 3
P1(I) = A(I) + F1 * (B(I) - A(I))
20 CONTINUE
30 CONTINUE
CALL LMRUN(VB, NV, P0, P1, IV0, IDSH)
RETURN
END
C
SUBROUTINE LMRUN(VB, NV, P0, P1, IV, IDSH)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), P0(3), P1(3)
INTEGER NV, IV, IDSH
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IV .EQ. 1) THEN
ISTYLE = 1
CALL MSEG(VB, NV, P0, P1)
ELSE IF (IDSH .EQ. 1) THEN
ISTYLE = 2
CALL MSEG(VB, NV, P0, P1)
ISTYLE = 1
END IF
C RESTOMOD END
RETURN
END
C
C LMOCC: 1 if the sight line to P passes through a placed solid
C other than IS (Cyrus-Beck against the face planes). Solids with
C the camera inside do not count (MPLACE).
INTEGER FUNCTION LMOCC(P, IS)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION P(3), T0, T1, DEN, T
INTEGER IS, K, J
C RESTOMOD BEGIN: Cyrus-Beck ray/convex test, published 1978
LMOCC = 0
DO 20 K = 1, NSOL
IF (K .EQ. IS) GO TO 20
IF (LACT(K) .EQ. 0 .OR. LINS(K) .EQ. 1) GO TO 20
T0 = 0.0D0
T1 = 0.999D0
DO 10 J = 1, NLF(K)
DEN = LWN(1,J,K)*P(1) + LWN(2,J,K)*P(2) + LWN(3,J,K)*P(3)
IF (DEN .EQ. 0.0D0) THEN
IF (LWD(J,K) .LT. 0.0D0) GO TO 20
ELSE
T = LWD(J,K) / DEN
IF (DEN .GT. 0.0D0) THEN
IF (T .LT. T1) T1 = T
ELSE
IF (T .GT. T0) T0 = T
END IF
END IF
IF (T0 .GE. T1) GO TO 20
10 CONTINUE
LMOCC = 1
RETURN
20 CONTINUE
C RESTOMOD END
RETURN
END
C
C-----------------------------------------------------------------------
C STKPL: place the docked CSM and LM. AT: the CSM's body axes in
C EQ; P (km, camera relative): where the centre of the CM's base
C is. The LM (model KL: KLMD gear down, KLMS stowed) faces it,
C its X axis against the CSM's and its Z axis along the CSM's (a
C half turn about Z; the roll between them is ours), tunnel top
C to tunnel top on the CSM's axis. The two tunnels' tops meet at
C the CSM tunnel's top, 10 ft 7 in above the CM's base (CSMBLD).
C-----------------------------------------------------------------------
SUBROUTINE STKPL(AT, P, KL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION AT(3,3), P(3), AL(3,3), PL(3), BO(3), Z(3)
INTEGER KL, I
CALL SETV(Z, 0.0D0, 0.0D0, 0.0D0)
CALL MPLACE(KCSM, AT, P, Z)
DO 10 I = 1, 3
AL(I,1) = -AT(I,1)
AL(I,2) = -AT(I,2)
AL(I,3) = AT(I,3)
PL(I) = P(I) + (10.0D0 + 7.0D0 / 12.0D0) * 0.3048D-3 * AT(I,1)
10 CONTINUE
CALL SETV(BO, 4.35D0, 0.0D0, -0.6D0)
CALL MPLACE(KL, AL, PL, BO)
RETURN
END
src/lvlab.f
C=======================================================================
C
C V I E W - 1 1 0 8 VEHICLE LABELS AND MARKERS
C
C Part of layer 6: called by MDRALL after the placed models are
C drawn. One relocatable element of the kernel; see vdrive.f for
C the list.
C
C A modern addition throughout (ours). VIEW drew "the vehicle
C outlines of the CSM, LM, and the S-IVB" at their apparent size
C (TN D-6853, printed p. 12); we have no source that it lettered
C them or marked the ones too small to see. Drawn only with a
C label level set (in_lablv 1-3), so in_lablv 0 keeps the picture
C it had.
C Labels (LB kind 8; id 1 CM, 2 SM, 3 LM, 4 S-IVB, 5 CSM) beside
C each placed model: off the model's projected X axis by the
C model's projected half-width across that axis plus a gap,
C so the name clears the model's lines. The CSM gets CM and SM
C labels beside its two modules, on the same side, or one CSM
C label when those two would touch.
C Markers: a placed model whose picture spans less than 0.2
C percent of the field (under about one plot pixel) gets the
C small boxed X of the Moon view's landing site (DMOON6, BOXX)
C at its centre instead, with its label beside the box. So do
C vehicles known only by their state: the CSM in scenes 5 and
C 6 (VSTATE, replay or tape), and the LM in its modelled
C descent (LMDESC, the last 600 s before touchdown) where the
C camera does not ride it. The LM has no state of its own
C anywhere else: the scenario's legs and the tape carry the CSM
C only. Docked, it rides with the CSM's mark; from undocking
C to the modelled descent, and after touchdown (where the Moon
C view's landing site mark stands), it is not marked.
C A label that would leave the frame or meet a name already
C lettered (the other layers' labels, VLSEED, or a vehicle's)
C tries the other side of its vehicle, then is dropped.
C=======================================================================
C
C-----------------------------------------------------------------------
C VLABEL: label the placed models and mark the vehicles too small
C to draw. Order CSM, LM, S-IVB: the first placed label wins.
C-----------------------------------------------------------------------
SUBROUTINE VLABEL(GET, VB, NV, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), LB(4,MAXL)
INTEGER NV, NL
DOUBLE PRECISION XA(64), XB(64), YA(64)
DOUBLE PRECISION YB(64), XM(8), YM(8), H, TINY, C(3), P(3), Q(3)
DOUBLE PRECISION R(3), V(3), XMN, XMX, YMN, YMX, X1, Y1, X2, Y2
DOUBLE PRECISION SALT, SEYE, FX(3), FY(3), FZ(3)
INTEGER KORD(4), IDM(4), NCH(5), M, K, N, J, NP, NP0, NM, IOK
INTEGER ISD
INTEGER LMKNOW
DATA KORD / KCSM, KLMD, KLMS, KSIV /
DATA IDM / 5, 3, 3, 4 /
DATA NCH / 2, 2, 2, 5, 3 /
NM = 0
C The name height TXALL letters at, and 0.2 percent of the field.
H = 0.028D0 * BOXH
TINY = 0.004D0 * BOXH
C The names already placed by the other layers are kept clear of.
CALL VLSEED(LB, NL, H, XA, XB, YA, YB, NP)
DO 50 M = 1, 4
K = KORD(M)
IF (MDON(K) .EQ. 0) GO TO 50
CALL VLPTS(K, VLPX, VLPY, N)
IF (N .EQ. 0) GO TO 50
XMN = VLPX(1)
XMX = VLPX(1)
YMN = VLPY(1)
YMX = VLPY(1)
DO 10 J = 2, N
XMN = DMIN1(XMN, VLPX(J))
XMX = DMAX1(XMX, VLPX(J))
YMN = DMIN1(YMN, VLPY(J))
YMX = DMAX1(YMX, VLPY(J))
10 CONTINUE
CALL VLCEN(K, MDS1(K), MDS2(K), C)
IF (DMAX1(XMX - XMN, YMX - YMN) .GE. TINY) GO TO 30
CALL VLMARK(VB, NV, LB, NL, C, IDM(M), NCH(IDM(M)), 0, H,
& XA, XB, YA, YB, NP, XM, YM, NM)
GO TO 50
C The CSM: the CM (its first two solids, cone and tunnel), then
C the SM (the rest) on the same side with the CM's box counted.
30 IF (K .NE. KCSM) GO TO 40
CALL VLCEN(K, MDS1(K), MDS1(K) + 1, P)
ISD = 0
CALL VLSIDE(K, P, VLPX, VLPY, N, 2, H, XA, XB, YA, YB, NP,
& X1, Y1, ISD, IOK)
IF (IOK .EQ. 0) GO TO 40
NP0 = NP
CALL VLADD(X1, Y1, 1.4D0 * H, H, XA, XB, YA, YB, NP)
CALL VLCEN(K, MDS1(K) + 2, MDS2(K), Q)
CALL VLSIDE(K, Q, VLPX, VLPY, N, 2, H, XA, XB, YA, YB, NP,
& X2, Y2, ISD, IOK)
NP = NP0
IF (IOK .EQ. 0) GO TO 40
CALL VLPUT(LB, NL, X1, Y1, H, 1, 2, XA, XB, YA, YB, NP)
CALL VLPUT(LB, NL, X2, Y2, H, 2, 2, XA, XB, YA, YB, NP)
GO TO 50
40 ISD = 0
CALL VLSIDE(K, C, VLPX, VLPY, N, NCH(IDM(M)), H, XA, XB, YA, YB,
& NP, X1, Y1, ISD, IOK)
IF (IOK .EQ. 1) CALL VLPUT(LB, NL, X1, Y1, H, IDM(M),
& NCH(IDM(M)), XA, XB, YA, YB, NP)
50 CONTINUE
C
C The CSM from its state, where it is not placed and the camera
C does not ride it (as TGTPOS): scenes 5 and 6.
IF (MDON(KCSM) .EQ. 1) GO TO 60
IF (ISCN .NE. 5 .AND. ISCN .NE. 6) GO TO 60
CALL VSTATE(GET, 2, R, V)
DO 55 J = 1, 3
P(J) = MPOS(J) + R(J)
55 CONTINUE
CALL VLMARK(VB, NV, LB, NL, P, 5, 3, 1, H, XA, XB, YA, YB, NP,
& XM, YM, NM)
C
C The LM in its modelled descent (LMDESC gives the commander's eye;
C LMDESC also sets LMALT and LMEYE for scene 5, kept here).
60 IF (LMKNOW(GET) .EQ. 0) RETURN
SALT = LMALT
SEYE = LMEYE
CALL LMDESC(GET, Q, FX, FY, FZ)
LMALT = SALT
LMEYE = SEYE
CALL MXV(MMF, Q, R)
DO 65 J = 1, 3
P(J) = MPOS(J) + R(J)
65 CONTINUE
CALL VLMARK(VB, NV, LB, NL, P, 3, 2, 1, H, XA, XB, YA, YB, NP,
& XM, YM, NM)
RETURN
END
C
C LMKNOW: 1 if the LM, not placed as a model and not carrying the
C camera (scene 5), has a state of its own at GET: the modelled
C descent, from 600 s before touchdown (LMDESC) to touchdown.
INTEGER FUNCTION LMKNOW(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET
LMKNOW = 0
IF (MDON(KLMD) .EQ. 1 .OR. MDON(KLMS) .EQ. 1) RETURN
IF (ISCN .EQ. 5 .OR. LUT0 .LE. 0.0D0) RETURN
IF (GET .LT. LUT0 - 600.0D0 .OR. GET .GE. LUT0) RETURN
LMKNOW = 1
RETURN
END
C
C VPRES: IVBIT, the vehicles in this frame's world, for hdr(21): 1
C the CSM, 2 the LM, 4 the S-IVB, each if placed as a model or
C known by its state for a marker (VLABEL); the vehicle the camera
C rides in a window or station view is not counted.
SUBROUTINE VPRES(GET)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET
INTEGER LMKNOW
IVBIT = 0
IF (MDON(KCSM) .EQ. 1 .OR. ISCN .EQ. 5 .OR. ISCN .EQ. 6)
& IVBIT = 1
IF (MDON(KLMD) .EQ. 1 .OR. MDON(KLMS) .EQ. 1 .OR.
& LMKNOW(GET) .EQ. 1) IVBIT = IVBIT + 2
IF (MDON(KSIV) .EQ. 1) IVBIT = IVBIT + 4
RETURN
END
C
C VLPTS: the plot points (PX, PY; N of them) of placed model K's
C solid vertices and free-line ends in front of the camera.
SUBROUTINE VLPTS(K, PX, PY, N)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K, N
DOUBLE PRECISION PX(2400), PY(2400), V(3), W(3), A(3)
INTEGER I, J, L, IS
N = 0
IF (MDS2(K) .LT. MDS1(K)) GO TO 30
DO 20 IS = MDS1(K), MDS2(K)
DO 10 J = 1, NLV(IS)
CALL VLPT1(LWV(1,J,IS), PX, PY, N)
10 CONTINUE
20 CONTINUE
30 IF (MDX2(K) .LT. MDX1(K)) RETURN
DO 50 J = MDX1(K), MDX2(K)
DO 45 L = 0, 1
DO 35 I = 1, 3
V(I) = (LXL(I + 3 * L, J) - MDBO(I,K)) * 1.0D-3
35 CONTINUE
CALL MXV(MDAT(1,1,K), V, W)
DO 40 I = 1, 3
A(I) = MDP(I,K) + W(I)
40 CONTINUE
CALL VLPT1(A, PX, PY, N)
45 CONTINUE
50 CONTINUE
RETURN
END
C
C VLPT1: add camera-relative point A's plot position, if it is in
C front of the camera and projected, to PX, PY.
SUBROUTINE VLPT1(A, PX, PY, N)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION A(3), PX(2400), PY(2400), X, Y
INTEGER N, IOK
IF (N .GE. 2400) RETURN
IF (A(1)*CB(1) + A(2)*CB(2) + A(3)*CB(3) .LE. 0.0D0) RETURN
CALL PROJ(A, X, Y, IOK)
IF (IOK .EQ. 0) RETURN
N = N + 1
PX(N) = X
PY(N) = Y
RETURN
END
C
C VLCEN: P (km, camera relative), the centre of the body-axis box
C around the vertices of placed model K's solids IS1..IS2.
SUBROUTINE VLCEN(K, IS1, IS2, P)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K, IS1, IS2
DOUBLE PRECISION P(3), BMN(3), BMX(3), V(3), W(3)
INTEGER I, J, IS
DO 10 I = 1, 3
BMN(I) = 1.0D30
BMX(I) = -1.0D30
10 CONTINUE
DO 30 IS = IS1, IS2
DO 25 J = 1, NLV(IS)
DO 20 I = 1, 3
BMN(I) = DMIN1(BMN(I), LMV(I,J,IS))
BMX(I) = DMAX1(BMX(I), LMV(I,J,IS))
20 CONTINUE
25 CONTINUE
30 CONTINUE
DO 40 I = 1, 3
V(I) = (0.5D0 * (BMN(I) + BMX(I)) - MDBO(I,K)) * 1.0D-3
40 CONTINUE
CALL MXV(MDAT(1,1,K), V, W)
DO 50 I = 1, 3
P(I) = MDP(I,K) + W(I)
50 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C VLSIDE: where to letter NC characters (height H) for the part of
C placed model K centred at P. The side is square to the model's
C projected X axis (right of it, or above it when the axis lies
C across the frame; right of the model when it is seen near end
C on); the text's near edge is the model's largest
C projected distance from that axis (over PX, PY) plus half a
C character height out. ISD 0 tries that side then the other
C and returns the one used (+1 or -1); ISD +-1 tries only that
C side. X0, Y0: the text's lower left. IOK = 1 if it fits in the
C frame clear of the labels placed so far (IVFREE).
C-----------------------------------------------------------------------
SUBROUTINE VLSIDE(K, P, PX, PY, N, NC, H, XA, XB, YA, YB, NP,
& X0, Y0, ISD, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K, N, NC, NP, ISD, IOK
DOUBLE PRECISION P(3), PX(2400), PY(2400), H, XA(64), XB(64)
DOUBLE PRECISION YA(64), YB(64), X0, Y0
DOUBLE PRECISION A(3), CX, CY, AX, AY, DX, DY, D, NX, NY, E, T
DOUBLE PRECISION WD, UX, UY, G, S
INTEGER I, J, L, L1, L2, IK, IVFREE
IOK = 0
IF (P(1)*CB(1) + P(2)*CB(2) + P(3)*CB(3) .LE. 0.0D0) RETURN
CALL PROJ(P, CX, CY, IK)
IF (IK .EQ. 0) RETURN
C The model's X axis in the plot, 1 m of it from P. Seen within
C 30 deg of end on, the axis says little: the label goes to the
C right (or left) of the whole model instead.
NX = 1.0D0
NY = 0.0D0
D = DSQRT(P(1) * P(1) + P(2) * P(2) + P(3) * P(3))
IF (DABS(P(1) * MDAT(1,1,K) + P(2) * MDAT(2,1,K)
& + P(3) * MDAT(3,1,K)) .GT. 0.866D0 * D) GO TO 15
DO 10 I = 1, 3
A(I) = P(I) + 1.0D-3 * MDAT(I,1,K)
10 CONTINUE
IF (A(1)*CB(1) + A(2)*CB(2) + A(3)*CB(3) .LE. 0.0D0) GO TO 15
CALL PROJ(A, AX, AY, IK)
DX = AX - CX
DY = AY - CY
D = DSQRT(DX * DX + DY * DY)
IF (IK .EQ. 0 .OR. D .LT. 1.0D-9) GO TO 15
NX = -DY / D
NY = DX / D
IF (NX .GT. 0.0D0 .OR. (NX .EQ. 0.0D0 .AND. NY .GT. 0.0D0))
& GO TO 15
NX = -NX
NY = -NY
C The model's half-width across the axis.
15 E = 0.0D0
DO 20 J = 1, N
T = DABS((PX(J) - CX) * NX + (PY(J) - CY) * NY)
IF (T .GT. E) E = T
20 CONTINUE
WD = 0.7D0 * H * DBLE(NC)
L1 = 1
L2 = 2
IF (ISD .EQ. 1) L2 = 1
IF (ISD .EQ. -1) L1 = 2
DO 40 L = L1, L2
S = DBLE(3 - 2 * L)
UX = S * NX
UY = S * NY
IF (DABS(UX) .LT. 0.5D0) GO TO 30
C Beside a steep axis: the text starts (or ends) G out.
G = E + 0.5D0 * H + 0.5D0 * H * DABS(UY)
X0 = CX + UX * G
IF (UX .LT. 0.0D0) X0 = X0 - WD
Y0 = CY + UY * G - 0.5D0 * H
GO TO 35
C Above or below a flat axis: the text centred G out.
30 G = E + 0.5D0 * H + 0.5D0 * WD * DABS(UX)
X0 = CX + UX * G - 0.5D0 * WD
Y0 = CY + UY * G
IF (UY .LT. 0.0D0) Y0 = Y0 - H
35 IF (IVFREE(X0, Y0, WD, H, XA, XB, YA, YB, NP) .EQ. 0) GO TO 40
IOK = 1
ISD = 3 - 2 * L
RETURN
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C VLMARK: mark a vehicle at P (km, camera relative) with a boxed X,
C and its label (id ID, NC characters, height H) at a corner of
C the box, the first of four that fits. Hidden behind the Earth
C or Moon; ISOL = 1 also hidden behind placed models and, from the
C LM window, below the sill (ISVIS mode 2), for a vehicle that is
C not itself placed. A second mark within a box of one already
C made (a docked stack) is left out.
C-----------------------------------------------------------------------
SUBROUTINE VLMARK(VB, NV, LB, NL, P, ID, NC, ISOL, H,
& XA, XB, YA, YB, NP, XM, YM, NM)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), LB(4,MAXL), P(3), H, XA(64), XB(64)
DOUBLE PRECISION YA(64), YB(64), XM(8), YM(8)
INTEGER NV, NL, ID, NC, ISOL, NP, NM
DOUBLE PRECISION X, Y, W, WD, X0, Y0, OCCL
INTEGER J, IOK, ISVIS, IVFREE
IOK = 1
IF (P(1)*CB(1) + P(2)*CB(2) + P(3)*CB(3) .LE. 0.0D0) RETURN
IF (OCCL(P, MPOS, RM) .GT. 0.0D0) RETURN
IF (OCCL(P, EPOS, RE) .GT. 0.0D0) RETURN
IVMODE = 2
IF (ISOL .EQ. 1) IOK = ISVIS(P)
IVMODE = 0
IF (ISOL .EQ. 1 .AND. IOK .EQ. 0) RETURN
CALL PROJ(P, X, Y, IOK)
IF (IOK .EQ. 0) RETURN
IF (DABS(X) .GT. BOXH .OR. DABS(Y) .GT. BOXH) RETURN
W = 0.012D0 * FOVH
DO 10 J = 1, NM
IF (DABS(X - XM(J)) .LT. W .AND. DABS(Y - YM(J)) .LT. W)
& RETURN
10 CONTINUE
CALL BOXX(VB, NV, X, Y, W)
IF (NM .GE. 8) GO TO 20
NM = NM + 1
XM(NM) = X
YM(NM) = Y
C The label at the box's upper right, upper left, lower right or
C lower left, the first that is free.
20 WD = 0.7D0 * H * DBLE(NC)
DO 25 J = 1, 4
X0 = X + W + 0.4D0 * H
IF (J .EQ. 2 .OR. J .EQ. 4) X0 = X - W - 0.4D0 * H - WD
Y0 = Y + W + 0.4D0 * H
IF (J .GE. 3) Y0 = Y - W - 0.4D0 * H - H
IF (IVFREE(X0, Y0, WD, H, XA, XB, YA, YB, NP) .EQ. 1) GO TO 30
25 CONTINUE
RETURN
30 CALL VLPUT(LB, NL, X0, Y0, H, ID, NC, XA, XB, YA, YB, NP)
RETURN
END
C
C IVFREE: 1 if a text box, lower left X0, Y0, WD wide and H high,
C lies in the frame (with its label point 0.4 H below and left, as
C TXALL letters it) and clear by 0.2 H of the NP boxes placed.
INTEGER FUNCTION IVFREE(X0, Y0, WD, H, XA, XB, YA, YB, NP)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION X0, Y0, WD, H, XA(64), XB(64), YA(64), YB(64)
DOUBLE PRECISION M
INTEGER NP, J
IVFREE = 0
IF (X0 - 0.4D0 * H .LT. -BOXH .OR. X0 + WD .GT. BOXH) RETURN
IF (Y0 - 0.4D0 * H .LT. -BOXH .OR. Y0 + H .GT. BOXH) RETURN
M = 0.2D0 * H
DO 10 J = 1, NP
IF (X0 - M .LT. XB(J) .AND. X0 + WD + M .GT. XA(J) .AND.
& Y0 - M .LT. YB(J) .AND. Y0 + H + M .GT. YA(J)) RETURN
10 CONTINUE
IVFREE = 1
RETURN
END
C
C VLSEED: the text boxes of the names the other layers have placed
C (LB kinds 1, 3-7, 9) that TXALL letters at this level, height H,
C as the first NP boxes. Crater names (kind 2) are the page's.
SUBROUTINE VLSEED(LB, NL, H, XA, XB, YA, YB, NP)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION LB(4,MAXL), H, XA(64), XB(64), YA(64), YB(64)
DOUBLE PRECISION X0, Y0
INTEGER NL, NP, I, J, K, ID, N
NP = 0
DO 50 I = 1, NL
K = NINT(LB(3,I))
ID = NINT(LB(4,I))
N = 0
IF ((K .EQ. 1 .OR. K .EQ. 6) .AND. ILEV .LT. 2) GO TO 50
IF (K .NE. 1 .OR. ID .LT. 1 .OR. ID .GT. NNAV) GO TO 15
DO 10 J = 1, 10
IF (NAVCH((ID - 1) * 10 + J) .NE. 0) N = N + 1
10 CONTINUE
15 IF (K .LT. 3 .OR. K .GT. 5) GO TO 25
DO 20 J = 1, 5
IF (BODCH((K - 3) * 5 + J) .NE. 0) N = N + 1
20 CONTINUE
25 IF (K .NE. 6 .OR. ID .LT. 1 .OR. ID .GT. NMARE) GO TO 35
DO 30 J = 1, 24
IF (MRCH((ID - 1) * 24 + J) .NE. 0) N = N + 1
30 CONTINUE
35 IF (K .EQ. 7) N = 22
IF (K .NE. 9) GO TO 45
DO 40 J = 1, 8
IF (PADCH(8 * (ISN - 1) + J) .NE. 0) N = N + 1
40 CONTINUE
45 IF (N .EQ. 0) GO TO 50
C Mare names are centred on the point, the others up and right.
X0 = LB(1,I) + 0.4D0 * H
Y0 = LB(2,I) + 0.4D0 * H
IF (K .EQ. 6) X0 = LB(1,I) - 0.35D0 * H * DBLE(N)
IF (K .EQ. 6) Y0 = LB(2,I) - 0.5D0 * H
CALL VLADD(X0, Y0, 0.7D0 * H * DBLE(N), H, XA, XB, YA, YB, NP)
50 CONTINUE
RETURN
END
C
C VLADD: count a text box (lower left X0, Y0, WD by H) as placed.
SUBROUTINE VLADD(X0, Y0, WD, H, XA, XB, YA, YB, NP)
DOUBLE PRECISION X0, Y0, WD, H, XA(64), XB(64), YA(64), YB(64)
INTEGER NP
IF (NP .GE. 64) RETURN
NP = NP + 1
XA(NP) = X0
XB(NP) = X0 + WD
YA(NP) = Y0
YB(NP) = Y0 + H
RETURN
END
C
C VLPUT: the vehicle label (LB kind 8, id ID) whose NC-character
C name TXALL letters with its lower left at X0, Y0: the label
C point is 0.4 H below and left of it.
SUBROUTINE VLPUT(LB, NL, X0, Y0, H, ID, NC, XA, XB, YA, YB, NP)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION LB(4,MAXL), X0, Y0, H, XA(64), XB(64), YA(64)
DOUBLE PRECISION YB(64)
INTEGER NL, ID, NC, NP
CALL LABEL(LB, NL, X0 - 0.4D0 * H, Y0 - 0.4D0 * H, 8, ID)
CALL VLADD(X0, Y0, 0.7D0 * H * DBLE(NC), H, XA, XB, YA, YB, NP)
RETURN
END
src/models.f
C=======================================================================
C
C V I E W - 1 1 0 8 SPACECRAFT MODELS
C
C Core element. The wireframe model library
C (data built once) that the vehicles layer draws. One
C relocatable element of the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C SPACECRAFT MODELS. A library of wireframe models in /CLM/,
C built once (MLIB). Each model is data: convex solids (prisms,
C MKPRS) and free lines (XLINE), in its own body frame, in metres.
C A free line with IS = 0 is a stand-alone line (legs, rims); with
C IS > 0 it is a mark on face LXF of solid IS (windows, target),
C hidden when that face is turned away. The model table /CMODI/
C gives each model's range of solids and lines.
C
C Per frame: MCLEAR, then MPLACE for each model in view (body axes,
C where a body point sits, camera relative), then MDRALL after the
C sky and bodies. Hidden parts: an edge between two faces turned
C away is hidden; any piece whose sight line passes through a
C placed solid is hidden (Cyrus-Beck, LMOCC); placed solids also
C hide stars, Sun, Earth and Moon (ISVIS, DSTARS, DSUN). Hidden
C pieces are dropped, as on the film, or drawn dashed (style 2)
C when IFLG bit 2 is set. TN D-6853 (printed p. 12): "Hidden-line
C models of the LM and the S-IVB can be produced".
C
C To add a model: a builder routine between MODBEG(K) and
C MODEND(K) in MLIB, a model number in viewcom.inc, and one MPLACE
C call in the scene's part of SCNMOD.
C=======================================================================
SUBROUTINE MLIB
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
NSOL = 0
NXL = 0
NMOD = 0
C LM, landing gear deployed (scene 4).
CALL MODBEG(KLMD)
CALL LMBODY
CALL LMGEAR(0)
CALL MODEND(KLMD)
C LM, landing gear stowed, with drogue and docking target, in its
C place on the S-IVB (scene 7).
CALL MODBEG(KLMS)
CALL LMBODY
CALL LMGEAR(1)
CALL LMDOCK
CALL MODEND(KLMS)
C S-IVB with the instrument unit and the stub of the SLA.
CALL MODBEG(KSIV)
CALL SIVBMD
CALL MODEND(KSIV)
C CSM, an outline (MDHL = 0).
CALL MODBEG(KCSM)
CALL CSMBLD
CALL MODEND(KCSM)
MDHL(KCSM) = 0
C CM cabin: the commander's window outlines about the design eye,
C an outline (station view only).
CALL MODBEG(KCMC)
CALL CMCAB
CALL MODEND(KCMC)
MDHL(KCMC) = 0
RETURN
END
C
C MODBEG, MODEND: open and close model K in the table.
SUBROUTINE MODBEG(K)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K
MDS1(K) = NSOL + 1
MDX1(K) = NXL + 1
MDHL(K) = 1
RETURN
END
C
SUBROUTINE MODEND(K)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K
MDS2(K) = NSOL
MDX2(K) = NXL
IF (K .GT. NMOD) NMOD = K
RETURN
END
C
C-----------------------------------------------------------------------
C LMBODY: the LM's stages, body axes X up, Y right, Z forward, the
C descent stage base at X = 0. Solids IS0+1 .. IS0+7; windows and
C hatch are marks on the cabin front (solid IS0+2, face 10 is its
C +Z cap).
C-----------------------------------------------------------------------
SUBROUTINE LMBODY
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION O(3), A1(3), A2(3), AN(3), P(2,8)
INTEGER IS0, IC
IS0 = NSOL
IC = IS0 + 2
C
C 1 descent stage: octagon 4.2 m across, X 0 .. 1.7.
CALL OCTAG(2.1D0, 2.1D0, 0.9D0, P)
CALL SETV(O, 0.0D0, 0.0D0, 0.0D0)
CALL SETV(A1, 0.0D0, 1.0D0, 0.0D0)
CALL SETV(A2, 0.0D0, 0.0D0, 1.0D0)
CALL SETV(AN, 1.0D0, 0.0D0, 0.0D0)
CALL MKPRS(8, P, O, A1, A2, AN, 1.7D0)
C 2 crew cabin, faceted, facing +Z.
CALL OCTAG(1.25D0, 1.0D0, 0.45D0, P)
CALL SETV(O, 2.9D0, 0.0D0, -0.2D0)
CALL SETV(A1, 0.0D0, 1.0D0, 0.0D0)
CALL SETV(A2, 1.0D0, 0.0D0, 0.0D0)
CALL SETV(AN, 0.0D0, 0.0D0, 1.0D0)
CALL MKPRS(8, P, O, A1, A2, AN, 1.35D0)
C 3 midsection behind the cabin.
CALL OCTAG(1.2D0, 1.05D0, 0.3D0, P)
CALL SETV(O, 2.95D0, 0.0D0, -1.5D0)
CALL MKPRS(8, P, O, A1, A2, AN, 1.3D0)
C 4 aft equipment bay.
CALL OCTAG(1.55D0, 0.5D0, 0.05D0, P)
CALL SETV(O, 3.1D0, 0.0D0, -2.2D0)
CALL MKPRS(8, P, O, A1, A2, AN, 0.7D0)
C 5 docking tunnel on top.
CALL OCTAG(0.5D0, 0.5D0, 0.2D0, P)
CALL SETV(O, 3.9D0, 0.0D0, -0.6D0)
CALL SETV(A1, 0.0D0, 1.0D0, 0.0D0)
CALL SETV(A2, 0.0D0, 0.0D0, 1.0D0)
CALL SETV(AN, 1.0D0, 0.0D0, 0.0D0)
CALL MKPRS(8, P, O, A1, A2, AN, 0.45D0)
C 6, 7 propellant tank bulges, left and right.
CALL OCTAG(0.6D0, 0.6D0, 0.25D0, P)
CALL SETV(O, 2.55D0, -1.2D0, -0.85D0)
CALL SETV(A1, 1.0D0, 0.0D0, 0.0D0)
CALL SETV(A2, 0.0D0, 0.0D0, 1.0D0)
CALL SETV(AN, 0.0D0, -1.0D0, 0.0D0)
CALL MKPRS(8, P, O, A1, A2, AN, 0.65D0)
CALL SETV(O, 2.55D0, 1.2D0, -0.85D0)
CALL SETV(AN, 0.0D0, 1.0D0, 0.0D0)
CALL MKPRS(8, P, O, A1, A2, AN, 0.65D0)
C
C Windows: two triangles on the cabin front.
CALL XLINE(-1.05D0, 3.55D0, 1.16D0, -0.35D0, 3.65D0, 1.16D0, IC)
CALL XLINE(-0.35D0, 3.65D0, 1.16D0, -0.55D0, 2.75D0, 1.16D0, IC)
CALL XLINE(-0.55D0, 2.75D0, 1.16D0, -1.05D0, 3.55D0, 1.16D0, IC)
CALL XLINE(1.05D0, 3.55D0, 1.16D0, 0.35D0, 3.65D0, 1.16D0, IC)
CALL XLINE(0.35D0, 3.65D0, 1.16D0, 0.55D0, 2.75D0, 1.16D0, IC)
CALL XLINE(0.55D0, 2.75D0, 1.16D0, 1.05D0, 3.55D0, 1.16D0, IC)
C Hatch on the cabin front.
CALL XLINE(-0.4D0, 2.0D0, 1.16D0, 0.4D0, 2.0D0, 1.16D0, IC)
CALL XLINE(0.4D0, 2.0D0, 1.16D0, 0.4D0, 2.6D0, 1.16D0, IC)
CALL XLINE(0.4D0, 2.6D0, 1.16D0, -0.4D0, 2.6D0, 1.16D0, IC)
CALL XLINE(-0.4D0, 2.6D0, 1.16D0, -0.4D0, 2.0D0, 1.16D0, IC)
RETURN
END
C
C-----------------------------------------------------------------------
C LMGEAR: landing gear as free lines. ISTOW = 0 deployed: on the
C diagonals a primary strut, two secondaries and a pad.
C ISTOW = 1 stowed. "In a retracted position until after the
C crew mans the LM, the landing gear struts are explosively
C extended" (Apollo 11 press kit, NASA release 69-83K, printed
C p. 103). How the folded gear lay is our guess: each primary
C strut runs up from its outrigger to a pad beside the ascent
C stage, the pad (37 in across, same page) square to the LM X
C axis; the secondaries are left out.
C-----------------------------------------------------------------------
SUBROUTINE LMGEAR(ISTOW)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER ISTOW
DOUBLE PRECISION SX, SZ, R0, R1, X0, X1, Q, C, S, C2, S2
INTEGER I, K
IF (ISTOW .EQ. 1) GO TO 30
DO 20 K = 0, 3
C = DCOS((45.0D0 + 90.0D0 * DBLE(K)) * DR)
S = DSIN((45.0D0 + 90.0D0 * DBLE(K)) * DR)
R0 = 2.33D0
R1 = 4.3D0
X0 = 1.5D0
X1 = -1.0D0
CALL XLINE(R0 * C, X0, R0 * S, R1 * C, X1, R1 * S, 0)
Q = 3.5D0
SX = Q * C
SZ = Q * S
CALL XLINE(SX, -0.45D0, SZ, 2.1D0 * DSIGN(1.0D0, C), 0.2D0,
& 1.2D0 * DSIGN(1.0D0, S), 0)
CALL XLINE(SX, -0.45D0, SZ, 1.2D0 * DSIGN(1.0D0, C), 0.2D0,
& 2.1D0 * DSIGN(1.0D0, S), 0)
DO 10 I = 0, 7
CALL XLINE(R1 * C + 0.45D0 * DCOS(DBLE(I) * PI / 4.0D0),
& X1 - 0.1D0, R1 * S + 0.45D0 * DSIN(DBLE(I) * PI / 4.0D0),
& R1 * C + 0.45D0 * DCOS(DBLE(I + 1) * PI / 4.0D0),
& X1 - 0.1D0,
& R1 * S + 0.45D0 * DSIN(DBLE(I + 1) * PI / 4.0D0), 0)
10 CONTINUE
20 CONTINUE
RETURN
30 DO 40 K = 0, 3
C = DCOS((45.0D0 + 90.0D0 * DBLE(K)) * DR)
S = DSIN((45.0D0 + 90.0D0 * DBLE(K)) * DR)
CALL XLINE(2.33D0 * C, 1.5D0, 2.33D0 * S,
& 2.3D0 * C, 2.2D0, 2.3D0 * S, 0)
DO 35 I = 0, 11
C2 = DCOS(DBLE(I) * PI / 6.0D0) * 0.47D0
S2 = DSIN(DBLE(I) * PI / 6.0D0) * 0.47D0
CALL XLINE(2.3D0 * C + C2, 2.2D0, 2.3D0 * S + S2,
& 2.3D0 * C + DCOS(DBLE(I + 1) * PI / 6.0D0) * 0.47D0, 2.2D0,
& 2.3D0 * S + DSIN(DBLE(I + 1) * PI / 6.0D0) * 0.47D0, 0)
35 CONTINUE
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C LMDOCK: docking drogue and CSM-active docking target, as marks
C on the LM model LMBODY has just built (tunnel = its solid 5,
C midsection = its solid 3).
C Drogue, a cone in the tunnel's +X cap (face 10): "a conical
C drogue mounted in the LM docking tunnel", which is 32 in across
C (press kit, printed p. 88 and p. 101). Depth: ours.
C Target for the CSM's crewman optical alignment sight (COAS): a
C disc on the midsection top (face 3) with a cross on a standoff
C above it, set off from the tunnel axis by as much as the COAS
C line of sight is from the CSM's docking axis (S7POSE), so it
C sits on the boresight. Size, standoff and offset: our guess.
C-----------------------------------------------------------------------
SUBROUTINE LMDOCK
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION C, S, C2, S2
INTEGER K, IT, IM
IT = NSOL - 2
IM = NSOL - 4
DO 70 K = 0, 11
C = DCOS(DBLE(K) * PI / 6.0D0)
S = DSIN(DBLE(K) * PI / 6.0D0)
C2 = DCOS(DBLE(K + 1) * PI / 6.0D0)
S2 = DSIN(DBLE(K + 1) * PI / 6.0D0)
CALL XLINE(0.4D0 * C, 4.35D0, -0.6D0 + 0.4D0 * S,
& 0.4D0 * C2, 4.35D0, -0.6D0 + 0.4D0 * S2, IT)
CALL XLINE(0.06D0 * C, 4.05D0, -0.6D0 + 0.06D0 * S,
& 0.06D0 * C2, 4.05D0, -0.6D0 + 0.06D0 * S2, IT)
IF (MOD(K, 3) .EQ. 0) CALL XLINE(0.4D0 * C, 4.35D0,
& -0.6D0 + 0.4D0 * S, 0.06D0 * C, 4.05D0, -0.6D0 + 0.06D0 * S,
& IT)
70 CONTINUE
DO 80 K = 0, 11
C = 0.18D0 * DCOS(DBLE(K) * PI / 6.0D0)
S = 0.18D0 * DSIN(DBLE(K) * PI / 6.0D0)
C2 = 0.18D0 * DCOS(DBLE(K + 1) * PI / 6.0D0)
S2 = 0.18D0 * DSIN(DBLE(K + 1) * PI / 6.0D0)
CALL XLINE(-0.72D0 + C, 4.0D0, -0.6D0 + S,
& -0.72D0 + C2, 4.0D0, -0.6D0 + S2, IM)
LXF(NXL) = 3
80 CONTINUE
CALL XLINE(-0.82D0, 4.45D0, -0.6D0, -0.62D0, 4.45D0, -0.6D0, 0)
CALL XLINE(-0.72D0, 4.45D0, -0.7D0, -0.72D0, 4.45D0, -0.5D0, 0)
RETURN
END
C
C-----------------------------------------------------------------------
C SIVBMD: S-IVB with the instrument unit on top, body axes X
C forward along the stage, origin at the centre of the top of the
C IU. One prism of 24 sides: both 21.7 ft across, 58.3 ft and
C 3 ft high (Apollo 11 press kit, printed p. 109).
C The SLA's fixed lower ring, left on the IU when the four upper
C panels were jettisoned at separation (AS-506 launch vehicle
C flight evaluation report, MPR-SAT-FE-69-9, p. xxiii; Apollo 11
C Flight Journal, 003:18:19). The SLA "is a truncated cone 28
C feet long tapering from 260 inches diameter at the base to 154
C inches at the forward end" (press kit, printed p. 88). The
C fixed panels are 7 ft high (secondary source: Wikipedia, "Apollo
C (spacecraft)"), so the ring's top is 233.5 in across. The
C jettisoned panels are not drawn; by the approach they had drifted
C off (our choice). Rim and four panel joints, the joints on the
C Y and Z axes (our guess).
C-----------------------------------------------------------------------
SUBROUTINE SIVBMD
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION O(3), A1(3), A2(3), AN(3), P24(2,24)
DOUBLE PRECISION RIU, XSL, RSL, Q, C, S, C2, S2
INTEGER K
RIU = 0.5D0 * 260.0D0 * 0.0254D0
DO 50 K = 1, 24
P24(1,K) = RIU * DCOS(DBLE(K) * PI / 12.0D0)
P24(2,K) = RIU * DSIN(DBLE(K) * PI / 12.0D0)
50 CONTINUE
Q = (58.3D0 + 3.0D0) * 0.3048D0
CALL SETV(O, -Q, 0.0D0, 0.0D0)
CALL SETV(A1, 0.0D0, 1.0D0, 0.0D0)
CALL SETV(A2, 0.0D0, 0.0D0, 1.0D0)
CALL SETV(AN, 1.0D0, 0.0D0, 0.0D0)
CALL MKPRS(24, P24, O, A1, A2, AN, Q)
XSL = 7.0D0 * 0.3048D0
RSL = 0.5D0 * (260.0D0 - 106.0D0 * 7.0D0 / 28.0D0) * 0.0254D0
DO 60 K = 0, 23
C = DCOS(DBLE(K) * PI / 12.0D0)
S = DSIN(DBLE(K) * PI / 12.0D0)
C2 = DCOS(DBLE(K + 1) * PI / 12.0D0)
S2 = DSIN(DBLE(K + 1) * PI / 12.0D0)
CALL XLINE(RSL * C, XSL, RSL * S, RSL * C2, XSL, RSL * S2, 0)
IF (MOD(K, 6) .EQ. 0) CALL XLINE(0.99D0 * RIU * C, 0.01D0,
& 0.99D0 * RIU * S, RSL * C, XSL, RSL * S, 0)
60 CONTINUE
RETURN
END
C
C XLINE: free line (IS=0) or mark on face LXF (10 unless set
C after the call: the cap of an 8-sided prism) of solid IS, given
C as Y, X, Z in the model's body metres.
SUBROUTINE XLINE(Y1, X1, Z1, Y2, X2, Z2, IS)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION Y1, X1, Z1, Y2, X2, Z2
INTEGER IS
IF (NXL .GE. MXL) RETURN
NXL = NXL + 1
LXL(1,NXL) = X1
LXL(2,NXL) = Y1
LXL(3,NXL) = Z1
LXL(4,NXL) = X2
LXL(5,NXL) = Y2
LXL(6,NXL) = Z2
LXS(NXL) = IS
LXF(NXL) = 10
RETURN
END
C
C OCTAG: chamfered rectangle, half sizes A by B, chamfer C, as 8
C points counter-clockwise (duplicates when C = 0 are harmless).
SUBROUTINE OCTAG(A, B, C, P)
DOUBLE PRECISION A, B, C, P(2,8)
P(1,1) = A
P(2,1) = -B + C
P(1,2) = A
P(2,2) = B - C
P(1,3) = A - C
P(2,3) = B
P(1,4) = -A + C
P(2,4) = B
P(1,5) = -A
P(2,5) = B - C
P(1,6) = -A
P(2,6) = -B + C
P(1,7) = -A + C
P(2,7) = -B
P(1,8) = A - C
P(2,8) = -B
RETURN
END
C
C MKPRS: prism over polygon P (NP points in axes A1, A2 about O),
C extruded H along AN. Faces 1..NP sides, NP+1 base, NP+2 cap.
SUBROUTINE MKPRS(NP, P, O, A1, A2, AN, H)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER NP
DOUBLE PRECISION P(2,NP), O(3), A1(3), A2(3), AN(3), H
DOUBLE PRECISION CEN(3), E1(3), E2(3), N(3), D, VDOT
INTEGER I, K, K2, IS, IA, IB, IC, NE
NSOL = NSOL + 1
IS = NSOL
NLV(IS) = 2 * NP
NLF(IS) = NP + 2
DO 20 K = 1, NP
DO 10 I = 1, 3
LMV(I,K,IS) = O(I) + P(1,K) * A1(I) + P(2,K) * A2(I)
LMV(I,K+NP,IS) = LMV(I,K,IS) + H * AN(I)
10 CONTINUE
20 CONTINUE
DO 25 I = 1, 3
CEN(I) = O(I) + 0.5D0 * H * AN(I)
25 CONTINUE
C Face planes, oriented away from the centre.
DO 40 K = 1, NP + 2
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (K .LE. NP) THEN
K2 = MOD(K, NP) + 1
IA = K
IB = K2
IC = K + NP
ELSE IF (K .EQ. NP + 1) THEN
IA = 1
IB = 3
IC = 6
ELSE
IA = 1 + NP
IB = 3 + NP
IC = 6 + NP
END IF
C RESTOMOD END
DO 30 I = 1, 3
E1(I) = LMV(I,IB,IS) - LMV(I,IA,IS)
E2(I) = LMV(I,IC,IS) - LMV(I,IA,IS)
30 CONTINUE
CALL VCRS(E1, E2, N)
CALL VUNIT(N)
D = VDOT(N, LMV(1,IA,IS))
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (D - VDOT(N, CEN) .LT. 0.0D0) THEN
DO 35 I = 1, 3
N(I) = -N(I)
35 CONTINUE
D = -D
END IF
C RESTOMOD END
DO 38 I = 1, 3
LMN(I,K,IS) = N(I)
38 CONTINUE
LMD(K,IS) = D
40 CONTINUE
C Edges: base, top, uprights.
NE = 0
DO 50 K = 1, NP
K2 = MOD(K, NP) + 1
NE = NE + 1
LME(1,NE,IS) = K
LME(2,NE,IS) = K2
LME(3,NE,IS) = K
LME(4,NE,IS) = NP + 1
NE = NE + 1
LME(1,NE,IS) = K + NP
LME(2,NE,IS) = K2 + NP
LME(3,NE,IS) = K
LME(4,NE,IS) = NP + 2
NE = NE + 1
LME(1,NE,IS) = K
LME(2,NE,IS) = K + NP
LME(3,NE,IS) = K
LME(4,NE,IS) = MOD(K + NP - 2, NP) + 1
50 CONTINUE
NLE(IS) = NE
RETURN
END
C
C-----------------------------------------------------------------------
C MKFRU: frustum over polygon P (NP points in axes A1, A2 about O),
C reaching H along AN with the top polygon scaled by SC; faces and
C edges numbered as MKPRS's (1..NP sides, NP+1 base, NP+2 top).
C-----------------------------------------------------------------------
SUBROUTINE MKFRU(NP, P, O, A1, A2, AN, H, SC)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER NP
DOUBLE PRECISION P(2,NP), O(3), A1(3), A2(3), AN(3), H, SC
DOUBLE PRECISION CEN(3), E1(3), E2(3), N(3), D, VDOT
INTEGER I, K, K2, IS, IA, IB, IC, NE
NSOL = NSOL + 1
IS = NSOL
NLV(IS) = 2 * NP
NLF(IS) = NP + 2
DO 20 K = 1, NP
DO 10 I = 1, 3
LMV(I,K,IS) = O(I) + P(1,K) * A1(I) + P(2,K) * A2(I)
LMV(I,K+NP,IS) = O(I) + H * AN(I)
& + SC * (P(1,K) * A1(I) + P(2,K) * A2(I))
10 CONTINUE
20 CONTINUE
DO 25 I = 1, 3
CEN(I) = O(I) + 0.5D0 * H * AN(I)
25 CONTINUE
DO 40 K = 1, NP + 2
K2 = MOD(K, NP) + 1
IA = K
IB = K2
IC = K + NP
IF (K .EQ. NP + 1) IA = 1
IF (K .EQ. NP + 1) IB = 3
IF (K .EQ. NP + 1) IC = 6
IF (K .EQ. NP + 2) IA = 1 + NP
IF (K .EQ. NP + 2) IB = 3 + NP
IF (K .EQ. NP + 2) IC = 6 + NP
DO 30 I = 1, 3
E1(I) = LMV(I,IB,IS) - LMV(I,IA,IS)
E2(I) = LMV(I,IC,IS) - LMV(I,IA,IS)
30 CONTINUE
CALL VCRS(E1, E2, N)
CALL VUNIT(N)
D = VDOT(N, LMV(1,IA,IS))
IF (D - VDOT(N, CEN) .GE. 0.0D0) GO TO 36
DO 35 I = 1, 3
N(I) = -N(I)
35 CONTINUE
D = -D
36 DO 38 I = 1, 3
LMN(I,K,IS) = N(I)
38 CONTINUE
LMD(K,IS) = D
40 CONTINUE
NE = 0
DO 50 K = 1, NP
K2 = MOD(K, NP) + 1
NE = NE + 1
LME(1,NE,IS) = K
LME(2,NE,IS) = K2
LME(3,NE,IS) = K
LME(4,NE,IS) = NP + 1
NE = NE + 1
LME(1,NE,IS) = K + NP
LME(2,NE,IS) = K2 + NP
LME(3,NE,IS) = K
LME(4,NE,IS) = NP + 2
NE = NE + 1
LME(1,NE,IS) = K
LME(2,NE,IS) = K + NP
LME(3,NE,IS) = K
LME(4,NE,IS) = MOD(K + NP - 2, NP) + 1
50 CONTINUE
NLE(IS) = NE
RETURN
END
C
C-----------------------------------------------------------------------
C CSMBLD: the command and service module. Body axes X along the
C stack toward the CM apex, Y and Z across (Apollo CSM convention:
C the RCS quads sit near +-Y and +-Z, below); origin at the centre
C of the CM's base, metres. Sources (CSM News Reference, North
C American Rockwell 1969, "NR"; Apollo Operations Handbook SM2A-
C 03-Block II-(1), 1969, "AOH"; Apollo 11 press kit, "PK"), where
C they disagree the one used is named:
C CM: "Height 10ft 7 in.", "Diameter 12ft 10in." (NR p. 39; PK
C p. 87 has 11 ft 5 in high, AOH p. 1-4 11 ft 1.5 in). The
C cone's shape is ours: a 16-sided frustum to 0.55 m radius at
C 2.25 m (about 32 deg half-angle), then a tunnel 0.45 m in
C radius to the 10 ft 7 in height (tunnel size ours).
C Fairing: "22 inches high" (NR p. 54; AOH p. 1-50 has 26 in),
C drawn as a short cylinder of the SM's diameter.
C SM: the cylinder "12 feet 11 inches long (high) and 12 feet 10
C inches in diameter" (AOH p. 1-50).
C SPS nozzle: an extension "protruding more than 9 feet below
C the aft bulkhead" (NR p. 58); we take 9 ft 8 in, the NR p. 3
C "22ft,7 in. excluding fairing" less the AOH's 12 ft 11 in
C (ours), and "an exit diameter of 7 feet 10-1/2 inches" (NR
C p. 162); the throat end (0.5 m radius) is ours.
C RCS quads: "four clusters of 90 degrees apart around the upper
C portion" (NR p. 58), "offset about 7 degrees from the Y and Z
C axes" (NR p. 148); each "eight feet long and nearly three feet
C wide" (NR p. 59), 0.3 m proud of the skin and centred 1.5 m
C below the SM's top (ours).
C High-gain antenna: on the aft bulkhead; "four 31-inch diameter
C parabolas" (NR p. 177); the boom "swings out at right angles
C to the spacecraft longitudinal axis, with the boom pointing 52
C degrees below the heads-up horizontal" (PK p. 90). Our
C reading: the boom leaves the SM's aft edge in the Y-Z plane,
C 52 deg from +Y toward -Z; its length (2.4 m) and the dishes'
C arrangement (2 by 2, square to the boom) are ours.
C-----------------------------------------------------------------------
SUBROUTINE CSMBLD
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION O(3), A1(3), A2(3), AN(3), P(2,24), Q(2,8)
DOUBLE PRECISION FT, RB, XF, XS, HS, XN, RN, C, S, C2, S2
DOUBLE PRECISION BR(3), BD(3), BU(3), T(3), CX, CY, HA, R
INTEGER K, J, I
FT = 0.3048D0
RB = 0.5D0 * (12.0D0 + 10.0D0 / 12.0D0) * FT
DO 10 K = 1, 16
P(1,K) = RB * DCOS(DBLE(K) * PI / 8.0D0)
P(2,K) = RB * DSIN(DBLE(K) * PI / 8.0D0)
10 CONTINUE
CALL SETV(A1, 0.0D0, 1.0D0, 0.0D0)
CALL SETV(A2, 0.0D0, 0.0D0, 1.0D0)
CALL SETV(AN, 1.0D0, 0.0D0, 0.0D0)
C CM cone and tunnel.
CALL SETV(O, 0.0D0, 0.0D0, 0.0D0)
CALL MKFRU(16, P, O, A1, A2, AN, 2.25D0, 0.55D0 / RB)
DO 12 K = 1, 16
P(1,K) = 0.45D0 * DCOS(DBLE(K) * PI / 8.0D0)
P(2,K) = 0.45D0 * DSIN(DBLE(K) * PI / 8.0D0)
12 CONTINUE
CALL SETV(O, 2.25D0, 0.0D0, 0.0D0)
CALL MKPRS(16, P, O, A1, A2, AN,
& (10.0D0 + 7.0D0 / 12.0D0) * FT - 2.25D0)
C SM and fairing, one cylinder of the SM's diameter up to the CM.
XF = 22.0D0 / 12.0D0 * FT
HS = (12.0D0 + 11.0D0 / 12.0D0) * FT
XS = -XF - HS
DO 14 K = 1, 16
P(1,K) = RB * DCOS(DBLE(K) * PI / 8.0D0)
P(2,K) = RB * DSIN(DBLE(K) * PI / 8.0D0)
14 CONTINUE
CALL SETV(O, XS, 0.0D0, 0.0D0)
CALL MKPRS(16, P, O, A1, A2, AN, HS + XF)
C The fairing's joint to the SM, a ring of free lines.
DO 16 K = 1, 16
CALL XLINE(P(1,K), -XF, P(2,K), P(1,MOD(K,16)+1), -XF,
& P(2,MOD(K,16)+1), 0)
16 CONTINUE
C SPS nozzle extension, from its throat end at the aft bulkhead.
XN = (9.0D0 + 8.0D0 / 12.0D0) * FT
RN = 0.5D0 * (7.0D0 + 10.5D0 / 12.0D0) * FT
DO 18 K = 1, 16
P(1,K) = 0.5D0 * DCOS(DBLE(K) * PI / 8.0D0)
P(2,K) = 0.5D0 * DSIN(DBLE(K) * PI / 8.0D0)
18 CONTINUE
CALL SETV(O, XS, 0.0D0, 0.0D0)
CALL SETV(AN, -1.0D0, 0.0D0, 0.0D0)
CALL MKFRU(16, P, O, A1, A2, AN, XN, RN / 0.5D0)
CALL SETV(AN, 1.0D0, 0.0D0, 0.0D0)
C RCS quad housings.
DO 30 J = 0, 3
C = DCOS((7.0D0 + 90.0D0 * DBLE(J)) * DR)
S = DSIN((7.0D0 + 90.0D0 * DBLE(J)) * DR)
CALL SETV(BR, 0.0D0, C, S)
CALL SETV(T, 0.0D0, -S, C)
CALL OCTAG(0.5D0 * 3.0D0 * FT, 0.15D0, 0.0D0, Q)
CALL SETV(O, -XF - 1.5D0 - 0.5D0 * 8.0D0 * FT,
& (RB + 0.15D0) * C, (RB + 0.15D0) * S)
CALL MKPRS(8, Q, O, T, BR, AN, 8.0D0 * FT)
30 CONTINUE
C High-gain antenna: boom and four dishes.
C = DCOS(-52.0D0 * DR)
S = DSIN(-52.0D0 * DR)
CALL SETV(BD, 0.0D0, C, S)
CALL SETV(BU, 1.0D0, 0.0D0, 0.0D0)
CALL VCRS(BD, BU, T)
CALL XLINE(RB * C, XS, RB * S, (RB + 2.4D0) * C, XS,
& (RB + 2.4D0) * S, 0)
HA = 0.5D0 * 31.0D0 * 0.0254D0
DO 50 I = 0, 3
CX = (DBLE(MOD(I, 2)) - 0.5D0) * 2.0D0 * HA
CY = (DBLE(I / 2) - 0.5D0) * 2.0D0 * HA
DO 40 K = 0, 11
C2 = DCOS(DBLE(K) * PI / 6.0D0) * HA
S2 = DSIN(DBLE(K) * PI / 6.0D0) * HA
R = RB + 2.4D0
CALL XLINE(R * C + (CX + C2) * T(2) + (CY + S2) * BU(2),
& XS + (CX + C2) * T(1) + (CY + S2) * BU(1),
& R * S + (CX + C2) * T(3) + (CY + S2) * BU(3),
& R * C + (CX + DCOS(DBLE(K + 1) * PI / 6.0D0) * HA) * T(2)
& + (CY + DSIN(DBLE(K + 1) * PI / 6.0D0) * HA) * BU(2),
& XS + (CX + DCOS(DBLE(K + 1) * PI / 6.0D0) * HA) * T(1)
& + (CY + DSIN(DBLE(K + 1) * PI / 6.0D0) * HA) * BU(1),
& R * S + (CX + DCOS(DBLE(K + 1) * PI / 6.0D0) * HA) * T(3)
& + (CY + DSIN(DBLE(K + 1) * PI / 6.0D0) * HA) * BU(3), 0)
40 CONTINUE
50 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C CMCAB: the CM cabin as the station view shows it: the left
C rendezvous window's outlines as VIEW drew them, in CSM body
C metres about the eye point CMEYE. MSC IN 69-FM-197, figure 9.0-3
C (PDF p. 263): "View as seen along CM X-axis during the PTC
C attitudes (X-axis in center of view)", gimbal angles 90, 0, 0,
C and "The CM left rendezvous window outline is shown for a zero
C roll attitude" (p. 17). The report draws two outlines; we read
C both off the plot by eye (to about 1 deg), in its plot degrees,
C where a point at plot radius R lies ATAN(R in radians) off the
C centre (as OVLPD). Our reading of the plot's axes: right is the
C CM's +Y, up its -Z (the rendezvous windows are on the -Z half,
C TN D-7439 p. 3), which puts the left window up and to the left.
C The x at the centre marks the CM X-axis, as in the report's CSM
C maneuver views ("the x also denotes the projection of the CM X-
C axis", MSC IN 69-FM-197 PDF p. 24, where "The CM left rendezvous
C window has been superimposed on these views"). The eye point is
C ours; the lines
C sit 0.5 m from it, so they are seen exactly along the outlines.
C Other windows: no outline we can read, so none drawn.
C-----------------------------------------------------------------------
SUBROUTINE CMCAB
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION WA(2,15), WB(2,13), E(3), A(3), B(3)
INTEGER K
DATA WA / -20.9D0, 11.5D0, -18.9D0, 16.9D0, -16.7D0, 22.6D0,
& -13.5D0, 28.4D0, -10.7D0, 33.4D0, -6.6D0, 38.6D0,
& -2.5D0, 38.8D0, -0.4D0, 34.5D0, 1.0D0, 30.5D0,
& 2.5D0, 25.5D0, 3.2D0, 20.5D0, -14.6D0, 4.3D0,
& -18.7D0, 8.3D0, -20.9D0, 10.4D0, -20.9D0, 11.5D0 /
DATA WB / -13.8D0, 10.4D0, -11.1D0, 16.9D0, -8.6D0, 21.9D0,
& -5.4D0, 26.9D0, -1.4D0, 32.7D0, 3.5D0, 38.4D0,
& 7.5D0, 38.8D0, 8.0D0, 34.5D0, 9.0D0, 34.2D0,
& 9.6D0, 29.1D0, 11.7D0, 19.0D0, -7.1D0, 4.3D0,
& -13.8D0, 10.4D0 /
CALL CMEYE(E)
DO 10 K = 1, 14
CALL CMDIR(E, WA(1,K), WA(2,K), A)
CALL CMDIR(E, WA(1,K+1), WA(2,K+1), B)
CALL XLINE(A(2), A(1), A(3), B(2), B(1), B(3), 0)
10 CONTINUE
DO 20 K = 1, 12
CALL CMDIR(E, WB(1,K), WB(2,K), A)
CALL CMDIR(E, WB(1,K+1), WB(2,K+1), B)
CALL XLINE(A(2), A(1), A(3), B(2), B(1), B(3), 0)
20 CONTINUE
CALL CMDIR(E, -0.7D0, -0.7D0, A)
CALL CMDIR(E, 0.7D0, 0.7D0, B)
CALL XLINE(A(2), A(1), A(3), B(2), B(1), B(3), 0)
CALL CMDIR(E, -0.7D0, 0.7D0, A)
CALL CMDIR(E, 0.7D0, -0.7D0, B)
CALL XLINE(A(2), A(1), A(3), B(2), B(1), B(3), 0)
RETURN
END
C
C CMEYE: the CM eye point, CSM body metres: 1.2 m above the CM's
C base, 0.5 m to -Y, 0.3 m to -Z (ours: no design-eye position
C found in our sources).
SUBROUTINE CMEYE(E)
DOUBLE PRECISION E(3)
E(1) = 1.2D0
E(2) = -0.5D0
E(3) = -0.3D0
RETURN
END
C
C CMDIR: the point 0.5 m from the eye E toward plot point (X, Y)
C of the station view (right +Y, up -Z, centre +X).
SUBROUTINE CMDIR(E, X, Y, P)
DOUBLE PRECISION E(3), X, Y, P(3), U(3), DR
DR = 3.141592653589793D0 / 180.0D0
U(1) = 1.0D0
U(2) = X * DR
U(3) = -Y * DR
CALL VUNIT(U)
P(1) = E(1) + 0.5D0 * U(1)
P(2) = E(2) + 0.5D0 * U(2)
P(3) = E(3) + 0.5D0 * U(3)
RETURN
END
src/pen.f
C=======================================================================
C
C V I E W - 1 1 0 8 PROJECTION AND PEN
C
C Core element. Directions to plot degrees,
C clipping, visibility, and vectors to the plot buffer. One
C relocatable element of the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C THE PEN. Points are camera-relative EQ vectors (km).
C PROJ direction to plot degrees
C UNPROJ plot degrees to a direction (reference or live camera)
C EMIT clip a plot-degree segment to the frame and store it
C SEG project and emit a 3-D segment
C PEN move (IP=0) or draw (IP=1) with the visibility test
C IVMODE; a segment crossing from seen to hidden is cut at
C the boundary by bisection.
C=======================================================================
C-----------------------------------------------------------------------
C PROJECTION. Radially symmetric about the boresight: a direction
C at angle T from the boresight and position angle P (from the
C right axis toward up) lands at radius RHO = K TAN(T/K), taken
C in degrees (times 180/PI, so RHO is T in degrees near the
C centre), at X = RHO COS P, Y = RHO SIN P. K = 1 (gnomonic, true
C perspective) up to a 100 deg field, rising linearly to K = 2
C (stereographic) at 170 deg (PK, set in VFRAME). The frame box
C half-width is RHO(FOV/2) (BOXH).
C Evidence, and it is our inference, not stated in either report:
C the film's descent frames (t28, t31, t35; FOV about 100) show
C the lunar horizon as a straight line at every height, and MSC
C IN 69-FM-197 (PDF p. 170, docking window, FOV 100) shows it
C straight across +-50 at Y = -30: a gnomonic plot draws every
C great circle straight. The same page's 170 deg front-window
C panel shows a fisheye dome, which a stereographic plot gives
C (conformal, circles stay circles); a gnomonic plot cannot reach
C 90 deg off axis at all. The film's evenly spaced ticks,
C labelled in degrees, fit a tangent-plane plot marked in degrees
C at the centre. (An earlier angle-angle mapping drew off-axis
C circles as rounded squares.) PTH returns T (rad).
C-----------------------------------------------------------------------
SUBROUTINE PROJ(D, X, Y, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION D(3), X, Y, A, B, C, S, T, R
INTEGER IOK
A = D(1)*CR(1) + D(2)*CR(2) + D(3)*CR(3)
B = D(1)*CU(1) + D(2)*CU(2) + D(3)*CU(3)
C = D(1)*CB(1) + D(2)*CB(2) + D(3)*CB(3)
S = DSQRT(A * A + B * B)
T = DATAN2(S, C)
PTH = T
IOK = 1
IF (T / DR .GT. THLIM) IOK = 0
IF (T / DR .GT. THLIM) T = THLIM * DR
R = PK * DTAN(T / PK) / DR
X = 0.0D0
Y = 0.0D0
IF (S .GT. 0.0D0) X = R * A / S
IF (S .GT. 0.0D0) Y = R * B / S
RETURN
END
C
C RHO: plot radius (deg) of a direction T (rad) off the boresight.
DOUBLE PRECISION FUNCTION RHO(T)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T
RHO = PK * DTAN(T / PK) / DR
RETURN
END
C
C UNPROJ: plot (X, Y) to a unit direction D. IREF=1: reference
C attitude with K = 1, for window overlays fixed to the vehicle
C (their plot coordinates were read off the film's 100 deg,
C gnomonic frames); IREF=0: the live camera and PK.
SUBROUTINE UNPROJ(X, Y, IREF, D)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION X, Y, D(3), A, B, C, R, T, Q
INTEGER IREF, I
R = DSQRT(X * X + Y * Y)
Q = PK
IF (IREF .EQ. 1) Q = 1.0D0
T = Q * DATAN(R * DR / Q)
C = DCOS(T)
A = 0.0D0
B = 0.0D0
IF (R .GT. 0.0D0) A = DSIN(T) * X / R
IF (R .GT. 0.0D0) B = DSIN(T) * Y / R
DO 10 I = 1, 3
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IREF .EQ. 1) THEN
D(I) = C * BREF(I) + A * RREF(I) + B * UREF(I)
ELSE
D(I) = C * CB(I) + A * CR(I) + B * CU(I)
END IF
C RESTOMOD END
10 CONTINUE
RETURN
END
C
C EMIT: Liang-Barsky clip to the frame, then store.
SUBROUTINE EMIT(VB, NV, X1, Y1, X2, Y2)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), X1, Y1, X2, Y2
INTEGER NV
DOUBLE PRECISION DX, DY, T0, T1, P(4), Q(4), R
INTEGER K
C RESTOMOD BEGIN: Liang-Barsky line clipping, published 1984
DX = X2 - X1
DY = Y2 - Y1
P(1) = -DX
Q(1) = X1 + BOXH
P(2) = DX
Q(2) = BOXH - X1
P(3) = -DY
Q(3) = Y1 + BOXH
P(4) = DY
Q(4) = BOXH - Y1
T0 = 0.0D0
T1 = 1.0D0
DO 10 K = 1, 4
IF (P(K) .EQ. 0.0D0) THEN
IF (Q(K) .LT. 0.0D0) RETURN
ELSE
R = Q(K) / P(K)
IF (P(K) .LT. 0.0D0) THEN
IF (R .GT. T1) RETURN
IF (R .GT. T0) T0 = R
ELSE
IF (R .LT. T0) RETURN
IF (R .LT. T1) T1 = R
END IF
END IF
10 CONTINUE
IF (NV .GE. MAXV) RETURN
NV = NV + 1
VB(1,NV) = X1 + T0 * DX
VB(2,NV) = Y1 + T0 * DY
VB(3,NV) = X1 + T1 * DX
VB(4,NV) = Y1 + T1 * DY
VB(5,NV) = DBLE(ISTYLE)
C RESTOMOD END
RETURN
END
C
C SEG: 3-D segment A-B (camera relative) to the frame. A segment
C running past the projection's limit (THLIM off the boresight) is
C cut at the limit by bisection.
SUBROUTINE SEG(VB, NV, A, B)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), A(3), B(3)
INTEGER NV
DOUBLE PRECISION X1, Y1, X2, Y2, P(3), Q(3), M(3), XM, YM
INTEGER K1, K2, KM, I, IT
CALL PROJ(A, X1, Y1, K1)
CALL PROJ(B, X2, Y2, K2)
IF (K1 .EQ. 0 .AND. K2 .EQ. 0) RETURN
IF (K1 .EQ. 1 .AND. K2 .EQ. 1) GO TO 50
C P inside the limit, Q outside.
DO 10 I = 1, 3
P(I) = A(I)
Q(I) = B(I)
IF (K1 .EQ. 0) P(I) = B(I)
IF (K1 .EQ. 0) Q(I) = A(I)
10 CONTINUE
DO 30 IT = 1, 20
DO 20 I = 1, 3
M(I) = 0.5D0 * (P(I) + Q(I))
20 CONTINUE
CALL PROJ(M, XM, YM, KM)
DO 25 I = 1, 3
IF (KM .EQ. 1) P(I) = M(I)
IF (KM .EQ. 0) Q(I) = M(I)
25 CONTINUE
30 CONTINUE
CALL PROJ(P, XM, YM, KM)
IF (K1 .EQ. 1) X2 = XM
IF (K1 .EQ. 1) Y2 = YM
IF (K1 .EQ. 0) X1 = XM
IF (K1 .EQ. 0) Y1 = YM
50 CALL EMIT(VB, NV, X1, Y1, X2, Y2)
RETURN
END
C
SUBROUTINE PEN(VB, NV, P, IP)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), P(3)
INTEGER NV, IP
DOUBLE PRECISION A(3), B(3), M(3)
INTEGER IV, ISVIS, I, IT
IV = ISVIS(P)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IP .EQ. 1) THEN
IF (IV .EQ. 1 .AND. IPV .EQ. 1) THEN
CALL SEG(VB, NV, PPX, P)
ELSE IF (IV .NE. IPV) THEN
C Seen end in A, hidden end in B; home in on the boundary.
DO 10 I = 1, 3
IF (IV .EQ. 1) THEN
A(I) = P(I)
B(I) = PPX(I)
ELSE
A(I) = PPX(I)
B(I) = P(I)
END IF
10 CONTINUE
DO 30 IT = 1, 10
DO 20 I = 1, 3
M(I) = 0.5D0 * (A(I) + B(I))
20 CONTINUE
IF (ISVIS(M) .EQ. 1) THEN
DO 22 I = 1, 3
A(I) = M(I)
22 CONTINUE
ELSE
DO 24 I = 1, 3
B(I) = M(I)
24 CONTINUE
END IF
30 CONTINUE
IF (IV .EQ. 1) THEN
CALL SEG(VB, NV, A, P)
ELSE
CALL SEG(VB, NV, PPX, A)
END IF
END IF
END IF
C RESTOMOD END
DO 40 I = 1, 3
PPX(I) = P(I)
40 CONTINUE
IPV = IV
RETURN
END
C
C-----------------------------------------------------------------------
C ISVIS: 1 if camera-relative point P is seen, by IVMODE:
C 0 always
C 1 on the Earth's surface: facing us, not behind the Moon
C 2 on the Earth's limb: not behind the Moon
C 3 on the Moon's surface: facing us, not behind the Earth
C 4 as 1, and on the night side
C 5 on the Moon's limb: not behind the Earth
C 6 on the Moon's surface, facing us and on the night side
C-----------------------------------------------------------------------
INTEGER FUNCTION ISVIS(P)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION P(3), N(3), OCCL, SILL
INTEGER I, LMOCC
ISVIS = 1
IF (IVMODE .EQ. 0) RETURN
C Placed spacecraft models hide what lies behind them.
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (NACT .GT. 0) THEN
IF (LMOCC(P, 0) .EQ. 1) THEN
ISVIS = 0
RETURN
END IF
END IF
C RESTOMOD END
C From the LM the window sill hides everything below it.
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF ((ISCN .EQ. 5 .AND. IVUSE .EQ. 0) .OR. IVUSE .EQ. 3) THEN
IF (SILL(P) .GT. 0.0D0) THEN
ISVIS = 0
RETURN
END IF
END IF
C RESTOMOD END
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IVMODE .EQ. 1 .OR. IVMODE .EQ. 4) THEN
DO 10 I = 1, 3
N(I) = P(I) - EPOS(I)
10 CONTINUE
IF (N(1)*P(1) + N(2)*P(2) + N(3)*P(3) .GE. 0.0D0) ISVIS = 0
IF (IVMODE .EQ. 4 .AND. ISVIS .EQ. 1) THEN
IF (N(1)*SUNU(1) + N(2)*SUNU(2) + N(3)*SUNU(3) .GT. 0.0D0)
& ISVIS = 0
END IF
IF (ISVIS .EQ. 1) THEN
IF (OCCL(P, MPOS, RM) .GT. 0.0D0) ISVIS = 0
END IF
ELSE IF (IVMODE .EQ. 2) THEN
IF (OCCL(P, MPOS, RM) .GT. 0.0D0) ISVIS = 0
ELSE IF (IVMODE .EQ. 3) THEN
DO 20 I = 1, 3
N(I) = P(I) - MPOS(I)
20 CONTINUE
IF (N(1)*P(1) + N(2)*P(2) + N(3)*P(3) .GE. 0.0D0) ISVIS = 0
IF (ISVIS .EQ. 1) THEN
IF (OCCL(P, EPOS, RE) .GT. 0.0D0) ISVIS = 0
END IF
ELSE IF (IVMODE .EQ. 5) THEN
IF (OCCL(P, EPOS, RE) .GT. 0.0D0) ISVIS = 0
ELSE IF (IVMODE .EQ. 6) THEN
DO 30 I = 1, 3
N(I) = P(I) - MPOS(I)
30 CONTINUE
IF (N(1)*P(1) + N(2)*P(2) + N(3)*P(3) .GE. 0.0D0) ISVIS = 0
IF (N(1)*SUNU(1) + N(2)*SUNU(2) + N(3)*SUNU(3) .GT. 0.0D0)
& ISVIS = 0
END IF
C RESTOMOD END
RETURN
END
C
C SILL: > 0 if P lies below the LM window sill, 35 deg below the
C centre of the reference (vehicle-fixed) frame.
DOUBLE PRECISION FUNCTION SILL(P)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION P(3), B, C
B = P(1)*UREF(1) + P(2)*UREF(2) + P(3)*UREF(3)
C = P(1)*BREF(1) + P(2)*BREF(2) + P(3)*BREF(3)
C Reference plot Y with K = 1 is TAN(T) SIN(P) = B / C (deg).
SILL = 1.0D0
IF (C .GT. 0.0D0) SILL = -35.0D0 - B / C / DR
IF (C .LE. 0.0D0 .AND. B .GE. 0.0D0) SILL = -1.0D0
RETURN
END
C
C OCCL: > 0 if the sight line from the camera to P passes inside
C the sphere of centre S (camera relative) and radius R.
DOUBLE PRECISION FUNCTION OCCL(P, S, R)
DOUBLE PRECISION P(3), S(3), R, T, PP, D1, D2, D3
PP = P(1)*P(1) + P(2)*P(2) + P(3)*P(3)
T = (P(1)*S(1) + P(2)*S(2) + P(3)*S(3)) / PP
IF (T .LT. 0.0D0) T = 0.0D0
IF (T .GT. 1.0D0) T = 1.0D0
D1 = T * P(1) - S(1)
D2 = T * P(2) - S(2)
D3 = T * P(3) - S(3)
OCCL = R * R - (D1 * D1 + D2 * D2 + D3 * D3)
RETURN
END
C
C RAYHIT: > 0 if the ray from the camera along unit U meets the
C sphere of centre S, radius R.
DOUBLE PRECISION FUNCTION RAYHIT(U, S, R)
DOUBLE PRECISION U(3), S(3), R, B
B = U(1)*S(1) + U(2)*S(2) + U(3)*S(3)
RAYHIT = -1.0D0
IF (B .LE. 0.0D0) RETURN
RAYHIT = R * R - (S(1)*S(1) + S(2)*S(2) + S(3)*S(3) - B * B)
RETURN
END
C
SUBROUTINE LABEL(LB, NL, X, Y, KIND, ID)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION LB(4,MAXL), X, Y
INTEGER NL, KIND, ID
IF (NL .GE. MAXL) RETURN
C With a label level set, the last 8 places are kept for the
C vehicle and pad labels (kinds 8, 9), which come after the sky's
C (ours: a full buffer of crater labels would crowd them out).
IF (ILABL .GE. 1 .AND. KIND .LT. 8 .AND. NL .GE. MAXL - 8) RETURN
IF (DABS(X) .GT. BOXH .OR. DABS(Y) .GT. BOXH) RETURN
NL = NL + 1
LB(1,NL) = X
LB(2,NL) = Y
LB(3,NL) = DBLE(KIND)
LB(4,NL) = DBLE(ID)
RETURN
END
C
C BOXX: a small boxed X of half-width W (plot deg) about (X, Y),
C the mark of the Moon view's landing site (DMOON6), also used for
C the launch pad (DPAD) and vehicle markers (VLMARK).
SUBROUTINE BOXX(VB, NV, X, Y, W)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), X, Y, W
INTEGER NV
CALL EMIT(VB, NV, X - W, Y - W, X + W, Y - W)
CALL EMIT(VB, NV, X + W, Y - W, X + W, Y + W)
CALL EMIT(VB, NV, X + W, Y + W, X - W, Y + W)
CALL EMIT(VB, NV, X - W, Y + W, X - W, Y - W)
CALL EMIT(VB, NV, X - W, Y - W, X + W, Y + W)
CALL EMIT(VB, NV, X - W, Y + W, X + W, Y - W)
RETURN
END
C
C-----------------------------------------------------------------------
C SHADE: night side of the sphere (centre S camera relative, radius
C R) as straight parallel lines (TN D-6853, printed p. 8, fig. 6).
C Each is cut from the sphere by a plane through the eye that
C contains the Sun's direction across the line of sight, so it
C draws as a straight line along the light; 15 span the disc.
C IVM is the visibility test for the night-side points.
C-----------------------------------------------------------------------
SUBROUTINE SHADE(VB, NV, S, R, IVM)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), S(3), R
INTEGER NV, IVM
DOUBLE PRECISION D, AE, U(3), Q(3), C(3), P(3), CF, LA, SF
INTEGER I, J
D = DSQRT(S(1)**2 + S(2)**2 + S(3)**2)
IF (D .LE. R) RETURN
AE = DASIN(R / D)
DO 10 I = 1, 3
U(I) = S(I) / D
10 CONTINUE
CF = SUNU(1)*U(1) + SUNU(2)*U(2) + SUNU(3)*U(3)
DO 20 I = 1, 3
Q(I) = SUNU(I) - CF * U(I)
20 CONTINUE
IF (Q(1)**2 + Q(2)**2 + Q(3)**2 .LT. 1.0D-12) RETURN
CALL VUNIT(Q)
CALL VCRS(U, Q, C)
IVMODE = IVM
DO 40 J = -7, 7
LA = DBLE(J) * AE / 7.5D0
DO 30 I = 1, 3
P(I) = C(I) * DCOS(LA) - U(I) * DSIN(LA)
30 CONTINUE
SF = D * DSIN(LA) / R
IF (DABS(SF) .LT. 1.0D0)
& CALL CIRCLE(VB, NV, S, R, P, DACOS(SF), 180)
40 CONTINUE
IVMODE = 0
RETURN
END
C
C-----------------------------------------------------------------------
C CIRCLE: small circle on the sphere (centre S camera relative,
C radius R), at angle G from the unit axis A, N points, drawn with
C the current IVMODE.
C-----------------------------------------------------------------------
SUBROUTINE CIRCLE(VB, NV, S, R, A, G, N)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), S(3), R, A(3), G
INTEGER NV, N
DOUBLE PRECISION E1(3), E2(3), P(3), T, CG, SG, CT, ST
INTEGER I, K
CALL PERP(A, E1, E2)
CG = DCOS(G)
SG = DSIN(G)
DO 20 K = 0, N
T = DBLE(K) * 2.0D0 * PI / DBLE(N)
CT = DCOS(T) * SG
ST = DSIN(T) * SG
DO 10 I = 1, 3
P(I) = S(I) + R * (CG * A(I) + CT * E1(I) + ST * E2(I))
10 CONTINUE
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (K .EQ. 0) THEN
CALL PEN(VB, NV, P, 0)
ELSE
CALL PEN(VB, NV, P, 1)
END IF
C RESTOMOD END
20 CONTINUE
RETURN
END
C
C PERP: unit vectors E1, E2 completing unit A to a right-handed set.
SUBROUTINE PERP(A, E1, E2)
DOUBLE PRECISION A(3), E1(3), E2(3), Z(3)
Z(1) = 0.0D0
Z(2) = 0.0D0
Z(3) = 1.0D0
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (DABS(A(3)) .GT. 0.9D0) THEN
Z(1) = 1.0D0
Z(3) = 0.0D0
END IF
C RESTOMOD END
CALL VCRS(Z, A, E1)
CALL VUNIT(E1)
CALL VCRS(A, E1, E2)
RETURN
END
C
C-----------------------------------------------------------------------
C MSEG: 3-D segment A-B (camera relative) to the frame, for model
C edges, which may pass beside or behind the camera or very close
C to it (a cabin seen from inside). Without recursion, a stack of
C pieces:
C both ends projected: one vector if the projected midpoint is
C within 0.1 percent of the box half-width of the chord (always
C so in the gnomonic plot, where lines stay straight), else
C split in two (a line bends in the stereographic plot, and
C can pass behind the camera between two seen ends);
C one end past the limit (THLIM, 90 deg times k): cut at the
C limit by bisection, keep the seen part;
C both ends past it: the limit cone is convex only up to 90 deg,
C so the piece can still cross the view: find its point
C nearest the boresight (the angle off it has one minimum along
C a line, so a ternary search) and split there if it is seen.
C Pieces are split at most 10 deep. A point at the camera itself
C projects to the centre (PROJ); nothing divides by the range.
C-----------------------------------------------------------------------
SUBROUTINE MSEG(VB, NV, A, B)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), A(3), B(3)
INTEGER NV
DOUBLE PRECISION SP(3,32), SQ(3,32), P(3), Q(3), M(3), T(3), W(3)
DOUBLE PRECISION XP, YP, XQ, YQ, XM, YM, DEV, TOL, T0, T1, TA, TB
DOUBLE PRECISION FA, FB, MSCOS
INTEGER SD(32), NS, D, KP, KQ, KM, I, IT
TOL = 1.0D-3 * BOXH
NS = 1
SD(1) = 0
DO 5 I = 1, 3
SP(I,1) = A(I)
SQ(I,1) = B(I)
5 CONTINUE
10 IF (NS .EQ. 0) RETURN
DO 12 I = 1, 3
P(I) = SP(I,NS)
Q(I) = SQ(I,NS)
12 CONTINUE
D = SD(NS)
NS = NS - 1
CALL PROJ(P, XP, YP, KP)
CALL PROJ(Q, XQ, YQ, KQ)
IF (KP .EQ. 0 .AND. KQ .EQ. 0) GO TO 40
IF (KP .EQ. 0 .OR. KQ .EQ. 0) GO TO 30
C Both ends seen.
DO 14 I = 1, 3
M(I) = 0.5D0 * (P(I) + Q(I))
14 CONTINUE
CALL PROJ(M, XM, YM, KM)
C Distance of the projected midpoint off the chord's line (the
C midpoint in space need not project to the chord's midpoint).
T0 = DSQRT((XQ - XP)**2 + (YQ - YP)**2)
DEV = DSQRT((XM - XP)**2 + (YM - YP)**2)
IF (T0 .GT. 1.0D-9) DEV = DABS((XM - XP) * (YQ - YP)
& - (YM - YP) * (XQ - XP)) / T0
IF (D .GE. 10) GO TO 20
IF (KM .EQ. 1 .AND. DEV .LE. TOL) GO TO 20
IF (NS .GE. 31) GO TO 20
CALL MSPUSH(SP, SQ, SD, NS, M, Q, D + 1)
CALL MSPUSH(SP, SQ, SD, NS, P, M, D + 1)
GO TO 10
20 CALL EMIT(VB, NV, XP, YP, XQ, YQ)
GO TO 10
C One end past the limit: bisect for the crossing (T seen, M not),
C keep the order A to B.
30 DO 32 I = 1, 3
T(I) = P(I)
M(I) = Q(I)
IF (KP .EQ. 0) T(I) = Q(I)
IF (KP .EQ. 0) M(I) = P(I)
32 CONTINUE
DO 36 IT = 1, 20
CALL MSMID(T, M, W, KM)
36 CONTINUE
IF (KP .EQ. 1) CALL MSPUSH(SP, SQ, SD, NS, P, T, D)
IF (KP .EQ. 0) CALL MSPUSH(SP, SQ, SD, NS, T, Q, D)
GO TO 10
C Both ends past the limit.
40 IF (D .GE. 10) GO TO 10
TA = 0.0D0
TB = 1.0D0
DO 44 IT = 1, 40
T0 = TA + (TB - TA) / 3.0D0
T1 = TB - (TB - TA) / 3.0D0
FA = MSCOS(P, Q, T0)
FB = MSCOS(P, Q, T1)
IF (FA .LT. FB) TA = T0
IF (FA .GE. FB) TB = T1
44 CONTINUE
T0 = 0.5D0 * (TA + TB)
DO 46 I = 1, 3
T(I) = P(I) + T0 * (Q(I) - P(I))
46 CONTINUE
CALL PROJ(T, XM, YM, KM)
IF (KM .EQ. 0 .OR. NS .GE. 31) GO TO 10
CALL MSPUSH(SP, SQ, SD, NS, T, Q, D + 1)
CALL MSPUSH(SP, SQ, SD, NS, P, T, D + 1)
GO TO 10
END
C
C MSPUSH: push the piece P-Q, depth D, on the MSEG stack.
SUBROUTINE MSPUSH(SP, SQ, SD, NS, P, Q, D)
DOUBLE PRECISION SP(3,32), SQ(3,32), P(3), Q(3)
INTEGER SD(32), NS, D, I
NS = NS + 1
DO 10 I = 1, 3
SP(I,NS) = P(I)
SQ(I,NS) = Q(I)
10 CONTINUE
SD(NS) = D
RETURN
END
C
C MSMID: one bisection step between seen T and unseen U; P is
C scratch. The midpoint replaces whichever end it matches.
SUBROUTINE MSMID(T, U, P, KM)
DOUBLE PRECISION T(3), U(3), P(3), X, Y
INTEGER KM, I
DO 10 I = 1, 3
P(I) = 0.5D0 * (T(I) + U(I))
10 CONTINUE
CALL PROJ(P, X, Y, KM)
DO 20 I = 1, 3
IF (KM .EQ. 1) T(I) = P(I)
IF (KM .EQ. 0) U(I) = P(I)
20 CONTINUE
RETURN
END
C
C MSCOS: cosine of the angle off the boresight of the point a
C fraction F from P to Q (-1 at the camera itself).
DOUBLE PRECISION FUNCTION MSCOS(P, Q, F)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION P(3), Q(3), F, X(3), R
INTEGER I
DO 10 I = 1, 3
X(I) = P(I) + F * (Q(I) - P(I))
10 CONTINUE
R = DSQRT(X(1) * X(1) + X(2) * X(2) + X(3) * X(3))
MSCOS = -1.0D0
IF (R .GT. 0.0D0) MSCOS = (X(1) * CB(1) + X(2) * CB(2)
& + X(3) * CB(3)) / R
RETURN
END
C
C OVLINE: straight line in reference plot degrees, carried to the
C live camera through its 3-D directions in 1 deg steps.
SUBROUTINE OVLINE(VB, NV, X1, Y1, X2, Y2)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION VB(5,MAXV), X1, Y1, X2, Y2
INTEGER NV
DOUBLE PRECISION D(3), F, L
INTEGER K, N
L = DSQRT((X2 - X1)**2 + (Y2 - Y1)**2)
N = 1 + INT(L)
DO 10 K = 0, N
F = DBLE(K) / DBLE(N)
CALL UNPROJ(X1 + F * (X2 - X1), Y1 + F * (Y2 - Y1), 1, D)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (K .EQ. 0) THEN
CALL PEN(VB, NV, D, 0)
ELSE
CALL PEN(VB, NV, D, 1)
END IF
C RESTOMOD END
10 CONTINUE
RETURN
END
src/sim.f
C=======================================================================
C
C V I E W - 1 1 0 8 THE ENGINE
C
C Core element. Flies the CSM through the current scenario: from
C its START state it integrates the motion, plays the SCORE (the
C BURN cards) as impulses, and at each REFERENCE row either resets
C the state to the sourced one (state vector updates on; our
C "delta correction") or only measures how far it has drifted
C (off). It writes the TAPE
C (tape.f) and knows nothing about drawing. One relocatable
C element of the kernel; see vdrive.f for the list.
C
C Period terms (docs/simulation.md). This plays the part of the
C RTACF integrator, "the Apollo Reference Mission Program", which
C wrote the trajectory ephemeris tape (Allday, TN D-6855, pp. 7-
C 8); the RTACF programs were taken over "without change to the
C basic logic and equations" and "run in a batch-processing mode"
C during missions (p. 7). A reset at a reference row plays a
C "CSM/LM STATE VECTOR UPDATE", uplinked to program P27 with verb
C 71 (Comanche 055, UPDATE_PROGRAM.agc). Fig. 2 of TN D-6855
C (p. 6) draws a "Vector transmit" line from the RTCC (IBM
C 360/75s, NASA-TM-X-64290, p. 113) to the RTACF, whose fig. 3
C computers are two Univac 1108s. The mapping is ours.
C VIEW's own integrator was Encke/Cowell (TN D-6853, p. 3).
C
C Physics (ours): Cowell's method, the CSM's geocentric position
C and velocity integrated directly, with the Earth, the Moon and
C the Sun as point masses (the Moon and Sun from ephem.f, the Moon
C through a half-hourly table, MOONQ; the Sun at 1 AU). Fourth-
C order Runge-Kutta; the step is 1/50 of the
C shorter of sqrt(r**3/mu) about the Earth and about the Moon,
C from 1 to 600 s, and it lands exactly on every burn and
C reference row. Burns are impulses at mid-burn. No Earth or
C Moon oblateness, no venting or attitude thrusting.
C
C=======================================================================
C
C-----------------------------------------------------------------------
C SIMRUN: run the engine over the current scenario and fill the
C tape. IFL bit 0: state vector updates on.
C-----------------------------------------------------------------------
SUBROUTINE SIMRUN(IFL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER IFL
DOUBLE PRECISION R(3), V(3), RR(3), VR(3), T, TEND, TN, H
DOUBLE PRECISION SIMDT
INTEGER K, J, I, IST, IDB(NBURN), IDR(NREF), ICOR
IF (ISN .EQ. 0) CALL SNSET(1)
CALL TPCLR
ISIMF = IFL
ITPSN = ISN
ICOR = MOD(IFL, 2)
DO 5 J = 1, NREF
RFOK(J) = 0
5 CONTINUE
IST = 0
DO 10 K = 1, NSTART
IF (STSN(K) .EQ. ISN) IST = K
10 CONTINUE
IF (IST .EQ. 0) RETURN
C Done flags: set for cards of other scenarios and unused slots.
DO 12 K = 1, NBURN
IDB(K) = 1
IF (K .LE. NBN) IDB(K) = 0
12 CONTINUE
DO 13 K = 1, NBN
IF (BNSN(K) .NE. ISN) IDB(K) = 1
13 CONTINUE
DO 14 K = 1, NREF
IDR(K) = 1
IF (K .LE. NRF) IDR(K) = 0
14 CONTINUE
DO 15 K = 1, NRF
IF (RFSN(K) .NE. ISN) IDR(K) = 1
15 CONTINUE
T = STP(3,IST)
TEND = STP(2,IST)
C The Moon at nodes every half hour over the run, for MOONQ.
MQDT = 1800.0D0
MQT0 = T - MQDT
NMQ = INT((TEND - MQT0) / MQDT) + 3
IF (NMQ .GT. MXMQ) NMQ = MXMQ
DO 16 K = 1, NMQ
CALL MOONV(MQT0 + DBLE(K - 1) * MQDT, MQ(1,K), MQ(4,K))
16 CONTINUE
CALL STATEV(STP(1,IST), STGC(IST), R, V)
CALL TPPUT(1, T, R, V)
C
C Step to the next event, or by the dynamic step.
20 IF (T .GE. TEND) RETURN
TN = TEND
DO 22 K = 1, NBURN
IF (IDB(K) .EQ. 0 .AND. BNT(K) .GT. T .AND. BNT(K) .LT. TN)
& TN = BNT(K)
22 CONTINUE
DO 24 K = 1, NREF
IF (IDR(K) .EQ. 0 .AND. RFP(3,K) .GT. T .AND. RFP(3,K) .LT. TN)
& TN = RFP(3,K)
24 CONTINUE
H = SIMDT(T, R)
IF (T + H .GE. TN) H = TN - T
CALL RK4(T, H, R, V)
T = T + H
IF (T .GE. TN) T = TN
CALL TPPUT(1, T, R, V)
IF (T .LT. TN) GO TO 20
C
C Reference rows at this time: measure, then correct if asked.
DO 40 J = 1, NREF
IF (IDR(J) .NE. 0 .OR. RFP(3,J) .NE. T) GO TO 40
IDR(J) = 1
CALL REFST(J, RR, VR)
RFERR(1,J) = DSQRT((R(1)-RR(1))**2 + (R(2)-RR(2))**2
& + (R(3)-RR(3))**2)
RFERR(2,J) = DSQRT((V(1)-VR(1))**2 + (V(2)-VR(2))**2
& + (V(3)-VR(3))**2) * 1.0D3 / 0.3048D0
RFOK(J) = 1
CALL TPMARK(T, 2, 1, J)
IF (ICOR .EQ. 0) GO TO 40
DO 30 I = 1, 3
R(I) = RR(I)
V(I) = VR(I)
30 CONTINUE
CALL TPPUT(1, T, R, V)
40 CONTINUE
C Burns at this time, as impulses.
DO 50 K = 1, NBURN
IF (IDB(K) .NE. 0 .OR. BNT(K) .NE. T) GO TO 50
IDB(K) = 1
CALL BURN(K, T, R, V)
CALL TPMARK(T, 1, 1, K)
CALL TPPUT(1, T, R, V)
50 CONTINUE
GO TO 20
END
C
C SIMDT: the integration step (s) at T, R: 1/50 of the local
C dynamical time about the nearer-acting body, 1 to 600 s.
DOUBLE PRECISION FUNCTION SIMDT(T, R)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T, R(3), PM(3), VM(3), RE3, RM3
CALL MOONQ(T, PM, VM)
RE3 = R(1)**2 + R(2)**2 + R(3)**2
RE3 = RE3 * DSQRT(RE3)
RM3 = (R(1)-PM(1))**2 + (R(2)-PM(2))**2 + (R(3)-PM(3))**2
RM3 = RM3 * DSQRT(RM3)
SIMDT = 0.02D0 * DMIN1(DSQRT(RE3 / GME), DSQRT(RM3 / GMM))
IF (SIMDT .LT. 1.0D0) SIMDT = 1.0D0
IF (SIMDT .GT. 600.0D0) SIMDT = 600.0D0
RETURN
END
C
C MOONQ: the Moon's geocentric position P and velocity V at T from
C the engine's node table (cubic Hermite; the series itself where
C T is off the table).
SUBROUTINE MOONQ(T, P, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T, P(3), V(3), S, H00, H10, H01, H11
DOUBLE PRECISION D00, D10, D01, D11
INTEGER K, I
K = INT((T - MQT0) / MQDT) + 1
IF (NMQ .LT. 2 .OR. K .LT. 1 .OR. K .GE. NMQ) GO TO 20
S = (T - MQT0) / MQDT - DBLE(K - 1)
H00 = (1.0D0 + 2.0D0 * S) * (1.0D0 - S)**2
H10 = S * (1.0D0 - S)**2
H01 = S * S * (3.0D0 - 2.0D0 * S)
H11 = S * S * (S - 1.0D0)
D00 = 6.0D0 * S * (S - 1.0D0)
D10 = (1.0D0 - S) * (1.0D0 - 3.0D0 * S)
D01 = -D00
D11 = S * (3.0D0 * S - 2.0D0)
DO 10 I = 1, 3
P(I) = H00 * MQ(I,K) + H10 * MQDT * MQ(I+3,K)
& + H01 * MQ(I,K+1) + H11 * MQDT * MQ(I+3,K+1)
V(I) = (D00 * MQ(I,K) + D10 * MQDT * MQ(I+3,K)
& + D01 * MQ(I,K+1) + D11 * MQDT * MQ(I+3,K+1)) / MQDT
10 CONTINUE
RETURN
20 CALL MOONV(T, P, V)
RETURN
END
C
C RK4: one Runge-Kutta step of H seconds from T for R, V.
SUBROUTINE RK4(T, H, R, V)
DOUBLE PRECISION T, H, R(3), V(3)
DOUBLE PRECISION A1(3), A2(3), A3(3), A4(3), R2(3), R3(3), R4(3)
DOUBLE PRECISION V2(3), V3(3), V4(3)
INTEGER I
CALL ACCEL(T, R, A1)
DO 10 I = 1, 3
R2(I) = R(I) + 0.5D0 * H * V(I)
V2(I) = V(I) + 0.5D0 * H * A1(I)
10 CONTINUE
CALL ACCEL(T + 0.5D0 * H, R2, A2)
DO 20 I = 1, 3
R3(I) = R(I) + 0.5D0 * H * V2(I)
V3(I) = V(I) + 0.5D0 * H * A2(I)
20 CONTINUE
CALL ACCEL(T + 0.5D0 * H, R3, A3)
DO 30 I = 1, 3
R4(I) = R(I) + H * V3(I)
V4(I) = V(I) + H * A3(I)
30 CONTINUE
CALL ACCEL(T + H, R4, A4)
DO 40 I = 1, 3
R(I) = R(I) + H / 6.0D0 * (V(I) + 2.0D0 * V2(I)
& + 2.0D0 * V3(I) + V4(I))
V(I) = V(I) + H / 6.0D0 * (A1(I) + 2.0D0 * A2(I)
& + 2.0D0 * A3(I) + A4(I))
40 CONTINUE
RETURN
END
C
C ACCEL: geocentric acceleration (km/s**2) at T, R: the Earth, and
C the Moon and the Sun as third bodies (direct minus indirect).
SUBROUTINE ACCEL(T, R, A)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T, R(3), A(3), PM(3), PS(3), D(3), R3, D3, P3
DOUBLE PRECISION GMS, AU, VM(3)
INTEGER I
GMS = 1.32712440018D11
AU = 1.495978707D8
CALL MOONQ(T, PM, VM)
CALL SUNG(T, PS)
R3 = R(1)**2 + R(2)**2 + R(3)**2
R3 = R3 * DSQRT(R3)
DO 10 I = 1, 3
A(I) = -GME * R(I) / R3
PS(I) = AU * PS(I)
D(I) = PM(I) - R(I)
10 CONTINUE
D3 = D(1)**2 + D(2)**2 + D(3)**2
D3 = D3 * DSQRT(D3)
P3 = PM(1)**2 + PM(2)**2 + PM(3)**2
P3 = P3 * DSQRT(P3)
DO 20 I = 1, 3
A(I) = A(I) + GMM * (D(I) / D3 - PM(I) / P3)
D(I) = PS(I) - R(I)
20 CONTINUE
D3 = D(1)**2 + D(2)**2 + D(3)**2
D3 = D3 * DSQRT(D3)
P3 = PS(1)**2 + PS(2)**2 + PS(3)**2
P3 = P3 * DSQRT(P3)
DO 30 I = 1, 3
A(I) = A(I) + GMS * (D(I) / D3 - PS(I) / P3)
30 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C BURN: apply burn K at T to R, V as an impulse: BNDV (ft/s) along
C BNP, BNR, BNN, the direction's components against the velocity,
C the in-plane radial and the orbit normal relative to the body
C BNBOD (1 Earth, 2 Moon).
C-----------------------------------------------------------------------
SUBROUTINE BURN(K, T, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER K, I
DOUBLE PRECISION T, R(3), V(3), RB(3), VB(3), PM(3), VM(3)
DOUBLE PRECISION X(3), Y(3), Z(3), DN, DV
DO 10 I = 1, 3
RB(I) = R(I)
VB(I) = V(I)
10 CONTINUE
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (BNBOD(K) .EQ. 2) THEN
CALL MOONV(T, PM, VM)
DO 20 I = 1, 3
RB(I) = R(I) - PM(I)
VB(I) = V(I) - VM(I)
20 CONTINUE
END IF
C RESTOMOD END
DO 30 I = 1, 3
X(I) = VB(I)
30 CONTINUE
CALL VUNIT(X)
CALL VCRS(RB, VB, Z)
CALL VUNIT(Z)
CALL VCRS(Z, X, Y)
DN = DSQRT(BNP(K)**2 + BNR(K)**2 + BNN(K)**2)
IF (DN .LE. 0.0D0) RETURN
DV = BNDV(K) * 0.3048D-3 / DN
DO 40 I = 1, 3
V(I) = V(I) + DV * (BNP(K) * X(I) + BNR(K) * Y(I)
& + BNN(K) * Z(I))
40 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C REFST: reference row J as a geocentric EQ state. About the Earth
C from its card as a CONIC leg's state (STATEV). About the Moon:
C its selenographic position at its altitude above the mean
C radius RM, its speed and flight-path angle, and the horizontal
C direction of the lunar leg's plane at that time (our choice),
C then the Moon's own state added.
C-----------------------------------------------------------------------
SUBROUTINE REFST(J, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER J, K, I, LEGAT
DOUBLE PRECISION R(3), V(3), S(3), U(3), M(3,3), HN(3), W(3)
DOUBLE PRECISION PM(3), VM(3), RA, SP, G
IF (RFBOD(J) .EQ. 2) GO TO 10
CALL STATEV(RFP(1,J), RFGC(J), R, V)
RETURN
10 CALL LLUNIT(RFP(4,J), RFP(5,J), S)
CALL MOONRT(RFP(3,J), M)
CALL MXV(M, S, U)
K = LEGAT(RFP(3,J), 2)
CALL VCRS(LGEL(1,K), LGEL(4,K), HN)
CALL VCRS(HN, U, W)
CALL VUNIT(W)
RA = RM + RFP(6,J) * 1.852D0
SP = RFP(7,J) * 0.3048D-3
G = RFP(8,J) * DR
CALL MOONV(RFP(3,J), PM, VM)
DO 20 I = 1, 3
R(I) = PM(I) + RA * U(I)
V(I) = VM(I) + SP * (DSIN(G) * U(I) + DCOS(G) * W(I))
20 CONTINUE
RETURN
END
src/tape.f
C=======================================================================
C
C V I E W - 1 1 0 8 THE TAPE
C
C Core element. Time-tagged vehicle states written by the engine
C (sim.f) and read by the state source (vsrc.f); it knows nothing
C of how they were made or how they are drawn. One relocatable
C element of the kernel; see vdrive.f for the list.
C
C Up to four vehicle channels (TN D-6853, printed p. 12: "As many
C as four vehicle trajectories can be integrated simultaneously");
C only channel 1, the CSM, is written now. A burn or a state
C vector update writes two samples at one time, before and after,
C and the reader never interpolates across such a pair. Reading
C is by cubic Hermite interpolation on the positions and
C velocities of the two samples around the time (ours).
C Period term: the "trajectory ephemeris tape", which "contained
C position and velocity-vector data and spacecraft-attitude
C information for one or two spacecraft" and was written by the
C RTACF integrator for programs "that contained no integrator"
C (Allday, TN D-6855, pp. 7-8). That VIEW read such a tape is our
C guess (TN D-6853 p. 3 names only the two parts).
C
C=======================================================================
C
C TPCLR: empty the tape.
SUBROUTINE TPCLR
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER C
DO 10 C = 1, MXCHN
NTP(C) = 0
10 CONTINUE
NMK = 0
RETURN
END
C
C TPPUT: append a sample of channel C at g.e.t. T, state R, V.
C A full channel keeps its last sample slot for the latest state.
SUBROUTINE TPPUT(C, T, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER C, I, K
DOUBLE PRECISION T, R(3), V(3)
K = NTP(C) + 1
IF (K .GT. MXSAM) K = MXSAM
NTP(C) = K
TPT(K,C) = T
DO 10 I = 1, 3
TPS(I,K,C) = R(I)
TPS(I+3,K,C) = V(I)
10 CONTINUE
RETURN
END
C
C TPMARK: append an event mark (kind KIND, channel C, index IX).
SUBROUTINE TPMARK(T, KIND, C, IX)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T
INTEGER KIND, C, IX
IF (NMK .GE. MXMRK) RETURN
NMK = NMK + 1
TMKT(NMK) = T
TMKK(NMK) = KIND
TMKC(NMK) = C
TMKI(NMK) = IX
RETURN
END
C
C-----------------------------------------------------------------------
C TPGET: state of channel C at g.e.t. T. IOK = 0 if T is outside
C the tape. At a time with two samples (a burn or an update)
C the later one is read.
C-----------------------------------------------------------------------
SUBROUTINE TPGET(C, T, R, V, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER C, IOK, LO, HI, MI, I
DOUBLE PRECISION T, R(3), V(3), H, S, H00, H10, H01, H11
DOUBLE PRECISION D00, D10, D01, D11
IOK = 0
IF (NTP(C) .LT. 2) RETURN
IF (T .LT. TPT(1,C) .OR. T .GT. TPT(NTP(C),C)) RETURN
IOK = 1
C The last sample LO with TPT(LO) <= T, by bisection.
LO = 1
HI = NTP(C)
10 IF (HI - LO .LE. 1) GO TO 20
MI = (LO + HI) / 2
IF (TPT(MI,C) .LE. T) LO = MI
IF (TPT(MI,C) .GT. T) HI = MI
GO TO 10
20 IF (TPT(HI,C) .LE. T) LO = HI
IF (LO .EQ. NTP(C)) LO = NTP(C) - 1
HI = LO + 1
H = TPT(HI,C) - TPT(LO,C)
C A zero-length step is a burn or an update: read its later side.
IF (H .GT. 0.0D0) GO TO 30
DO 25 I = 1, 3
R(I) = TPS(I,HI,C)
V(I) = TPS(I+3,HI,C)
25 CONTINUE
RETURN
30 S = (T - TPT(LO,C)) / H
IF (S .GT. 1.0D0) S = 1.0D0
C Cubic Hermite basis and its derivative (per unit S).
H00 = (1.0D0 + 2.0D0 * S) * (1.0D0 - S)**2
H10 = S * (1.0D0 - S)**2
H01 = S * S * (3.0D0 - 2.0D0 * S)
H11 = S * S * (S - 1.0D0)
D00 = 6.0D0 * S * (S - 1.0D0)
D10 = (1.0D0 - S) * (1.0D0 - 3.0D0 * S)
D01 = -D00
D11 = S * (3.0D0 * S - 2.0D0)
DO 40 I = 1, 3
R(I) = H00 * TPS(I,LO,C) + H10 * H * TPS(I+3,LO,C)
& + H01 * TPS(I,HI,C) + H11 * H * TPS(I+3,HI,C)
V(I) = (D00 * TPS(I,LO,C) + D10 * H * TPS(I+3,LO,C)
& + D01 * TPS(I,HI,C) + D11 * H * TPS(I+3,HI,C)) / H
40 CONTINUE
RETURN
END
src/traj.f
C=======================================================================
C
C V I E W - 1 1 0 8 TRAJECTORY LEGS
C
C Core element. Where the spacecraft is at a GET. One relocatable
C element of the kernel; see vdrive.f for the list.
C
C The scenario (data/scenarios, BLOCK DATA /CSCEN/: one mission's
C data, a run deck in 1969 terms) gives the
C epoch, the trajectory legs and the events. VIEW itself
C integrated trajectories from the state vectors of the operational
C trajectory document (TN D-6853, p. 3, 12); we put each leg on
C one simple model instead, fixed by sourced states:
C CIRC Earth circular orbit through a state (parking orbit);
C CONIC Earth-centred Kepler conic from a state, no lunar
C gravity (translunar and transearth coast);
C LUNAR circle about the Moon through two states, its plane
C and mean motion from those states (lunar orbit).
C A g.e.t. outside every leg takes the nearest leg.
C LM descent: P64 approach from 7200 ft altitude and 25600 ft
C range at GET 102:41:30 down to the site at touchdown.
C
C=======================================================================
C
C-----------------------------------------------------------------------
C SNSET: make scenario IM current: its epoch TJD0 (the ephemeris
C keys on it), the elements of each of its legs (LGEL), and the
C event times the scenes use.
C-----------------------------------------------------------------------
SUBROUTINE SNSET(IM)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER IM
DOUBLE PRECISION R(3), V(3), A(3), B(3), H(3), Y(3), M(3,3)
DOUBLE PRECISION S(3), E(3), EV, AX, RR, VV, CN, SN, NU, EA, AN
DOUBLE PRECISION EVGET, VDOT, VNRM, ANG
INTEGER K, I
ISN = IM
TJD0 = SNJD0(IM)
DO 90 K = 1, NLEG
IF (LGSN(K) .NE. IM) GO TO 90
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (LGTYP(K) .EQ. KCIRC) THEN
C Circle through the state's position, along its heading.
CALL STATEV(LGP(1,K), LGGC(K), R, V)
CALL VUNIT(R)
CALL VUNIT(V)
DO 10 I = 1, 3
LGEL(I,K) = R(I)
LGEL(I+3,K) = V(I)
10 CONTINUE
LGEL(9,K) = LGP(3,K)
LGEL(10,K) = RE + LGP(6,K) * 1.852D0
LGEL(11,K) = DSQRT(GME / LGEL(10,K)**3)
ELSE IF (LGTYP(K) .EQ. KCONIC) THEN
C Conic elements from the state: perigee unit P (1-3), Q 90
C deg ahead (4-6), eccentricity (7), semi-major axis (8),
C perigee time (9).
CALL STATEV(LGP(1,K), LGGC(K), R, V)
CALL VCRS(R, V, H)
RR = VNRM(R)
VV = VDOT(V, V)
CALL VCRS(V, H, E)
DO 20 I = 1, 3
E(I) = E(I) / GME - R(I) / RR
20 CONTINUE
EV = VNRM(E)
AX = 1.0D0 / (2.0D0 / RR - VV / GME)
CALL VUNIT(H)
CALL VUNIT(E)
CALL VCRS(H, E, Y)
CN = VDOT(R, E) / RR
SN = VDOT(R, Y) / RR
NU = DATAN2(SN, CN)
EA = 2.0D0 * DATAN(DSQRT((1.0D0 - EV) / (1.0D0 + EV))
& * DTAN(0.5D0 * NU))
AN = DSQRT(GME / AX**3)
DO 30 I = 1, 3
LGEL(I,K) = E(I)
LGEL(I+3,K) = Y(I)
30 CONTINUE
LGEL(7,K) = EV
LGEL(8,K) = AX
LGEL(9,K) = LGP(3,K) - (EA - EV * DSIN(EA)) / AN
ELSE
C Lunar circle through state A at T and state B at TB, both
C selenographic, carried to EQ by the Moon's orientation at
C their times. Sense: retrograde (LGGC = 0), the westward
C motion of MR Table 7-II's lunar rows, or prograde (1).
CALL LLUNIT(LGP(4,K), LGP(5,K), S)
CALL MOONRT(LGP(3,K), M)
CALL MXV(M, S, A)
CALL LLUNIT(LGP(11,K), LGP(12,K), S)
CALL MOONRT(LGP(10,K), M)
CALL MXV(M, S, B)
CALL VCRS(A, B, H)
CALL VUNIT(H)
DO 40 I = 1, 3
S(I) = M(I,3)
40 CONTINUE
IF ((VDOT(H, S) .GT. 0.0D0) .EQV. (LGGC(K) .EQ. 0)) THEN
DO 45 I = 1, 3
H(I) = -H(I)
45 CONTINUE
END IF
CALL VCRS(H, A, Y)
ANG = DATAN2(VDOT(B, Y), VDOT(B, A))
IF (ANG .LT. 0.0D0) ANG = ANG + 2.0D0 * PI
DO 50 I = 1, 3
LGEL(I,K) = A(I)
LGEL(I+3,K) = Y(I)
50 CONTINUE
LGEL(9,K) = LGP(3,K)
LGEL(10,K) = RM + LGP(6,K) * 1.852D0
LGEL(11,K) = (ANG + 2.0D0 * PI * DBLE(LGN(K)))
& / (LGP(10,K) - LGP(3,K))
END IF
C RESTOMOD END
90 CONTINUE
LUT0 = EVGET(KETD)
TETP = EVGET(KEEI)
RETURN
END
C
C EVGET: g.e.t. (s) of the current scenario's event of kind KIND,
C or -1 if it has none.
DOUBLE PRECISION FUNCTION EVGET(KIND)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER KIND, J
EVGET = -1.0D0
DO 10 J = 1, NEVT
IF (EVSN(J) .EQ. ISN .AND. EVKND(J) .EQ. KIND) EVGET = EVT(J)
10 CONTINUE
RETURN
END
C
C LEGAT: the current scenario's leg about the Earth (ICLS = 1: CIRC
C or CONIC) or the Moon (ICLS = 2: LUNAR) whose span holds GET,
C else the one whose span ends nearest it; 0 if there is none.
INTEGER FUNCTION LEGAT(GET, ICLS)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, D, DBEST
INTEGER ICLS, K, IC
LEGAT = 0
DBEST = 1.0D30
DO 10 K = 1, NLEG
IF (LGSN(K) .NE. ISN) GO TO 10
IC = 1
IF (LGTYP(K) .EQ. KLUNAR) IC = 2
IF (IC .NE. ICLS) GO TO 10
D = 0.0D0
IF (GET .LT. LGP(1,K)) D = LGP(1,K) - GET
IF (GET .GT. LGP(2,K)) D = GET - LGP(2,K)
IF (D .GE. DBEST) GO TO 10
DBEST = D
LEGAT = K
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C STATEV: the state in card fields P (the layout of LGP), geocentric
C EQ km and km/s. Latitude
C geodetic (MR Table 7-I, p. 7-8) or geocentric (IGC = 1, as in
C SP-4029's ascent table); altitude above the ellipsoid (ours:
C equatorial radius RE, flattening 1/298.257); longitude Earth
C fixed, turned by GMST at the state's time to the equator and
C equinox of date, then to J2000 by PRECM's transpose. Speed,
C flight-
C path angle and heading are space-fixed, against the geocentric
C horizontal (MR Table 7-I).
C-----------------------------------------------------------------------
SUBROUTINE STATEV(P, IGC, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER IGC
DOUBLE PRECISION P(NLGP), R(3), V(3)
DOUBLE PRECISION F, E2, FI, LA, H, SF, CF, XN, U(3), N(3), E(3)
DOUBLE PRECISION G, HD, SP, GMSTAT, PSI, RA, PM(3,3), W(3)
INTEGER I
F = 1.0D0 / 298.257D0
E2 = F * (2.0D0 - F)
FI = P(4) * DR
LA = P(5) * DR + GMSTAT(P(3))
H = P(6) * 1.852D0
SF = DSIN(FI)
CF = DCOS(FI)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IGC .EQ. 1) THEN
RA = RE * (1.0D0 - F * SF * SF) + H
R(1) = RA * CF * DCOS(LA)
R(2) = RA * CF * DSIN(LA)
R(3) = RA * SF
ELSE
XN = RE / DSQRT(1.0D0 - E2 * SF * SF)
R(1) = (XN + H) * CF * DCOS(LA)
R(2) = (XN + H) * CF * DSIN(LA)
R(3) = (XN * (1.0D0 - E2) + H) * SF
END IF
C RESTOMOD END
C Geocentric horizontal at R.
DO 10 I = 1, 3
U(I) = R(I)
10 CONTINUE
CALL VUNIT(U)
PSI = DASIN(U(3))
N(1) = -DSIN(PSI) * DCOS(LA)
N(2) = -DSIN(PSI) * DSIN(LA)
N(3) = DCOS(PSI)
E(1) = -DSIN(LA)
E(2) = DCOS(LA)
E(3) = 0.0D0
SP = P(7) * 0.3048D-3
G = P(8) * DR
HD = P(9) * DR
DO 20 I = 1, 3
V(I) = SP * (DSIN(G) * U(I) + DCOS(G) * (DCOS(HD) * N(I)
& + DSIN(HD) * E(I)))
20 CONTINUE
C Equator of date to J2000.
CALL PRECM((TJD0 + P(3) / 86400.0D0 - 2451545.0D0) / 36525.0D0,
& PM)
CALL MTXV(PM, R, W)
CALL MTXV(PM, V, U)
DO 30 I = 1, 3
R(I) = W(I)
V(I) = U(I)
30 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C ERTORB: spacecraft about the Earth, geocentric EQ km and km/s,
C on the current scenario's Earth leg for GET.
C-----------------------------------------------------------------------
SUBROUTINE ERTORB(GET, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, R(3), V(3), TH, C, S
INTEGER K, I, LEGAT
K = LEGAT(GET, 1)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (LGTYP(K) .EQ. KCONIC) THEN
CALL KEPLER(LGEL(1,K), LGEL(4,K), LGEL(7,K), LGEL(8,K),
& LGEL(9,K), GET, R, V)
ELSE
TH = LGEL(11,K) * (GET - LGEL(9,K))
C = DCOS(TH)
S = DSIN(TH)
DO 10 I = 1, 3
R(I) = LGEL(10,K) * (C * LGEL(I,K) + S * LGEL(I+3,K))
V(I) = LGEL(10,K) * LGEL(11,K)
& * (C * LGEL(I+3,K) - S * LGEL(I,K))
10 CONTINUE
END IF
C RESTOMOD END
RETURN
END
C
C-----------------------------------------------------------------------
C LUNORB: CSM relative to the Moon, EQ km and km/s, on the current
C scenario's lunar leg for GET.
C-----------------------------------------------------------------------
SUBROUTINE LUNORB(GET, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, R(3), V(3), TH, C, S
INTEGER K, I, LEGAT
K = LEGAT(GET, 2)
TH = LGEL(11,K) * (GET - LGEL(9,K))
C = DCOS(TH)
S = DSIN(TH)
DO 10 I = 1, 3
R(I) = LGEL(10,K) * (C * LGEL(I,K) + S * LGEL(I+3,K))
V(I) = LGEL(10,K) * LGEL(11,K)
& * (C * LGEL(I+3,K) - S * LGEL(I,K))
10 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C KEPLER: position and velocity at GET on the ellipse P, Q, E, A
C with perigee at time TP.
C-----------------------------------------------------------------------
SUBROUTINE KEPLER(P, Q, E, A, TP, GET, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION P(3), Q(3), E, A, TP, GET, R(3), V(3)
DOUBLE PRECISION AM, EA, F, RR, B, C, S, VF
INTEGER I, IT
AM = DSQRT(GME / A**3) * (GET - TP)
EA = AM
IF (E .GT. 0.8D0) EA = PI * DSIGN(1.0D0, AM)
IF (DABS(AM) .GT. PI) EA = AM
DO 10 IT = 1, 50
F = (EA - E * DSIN(EA) - AM) / (1.0D0 - E * DCOS(EA))
EA = EA - F
IF (DABS(F) .LT. 1.0D-12) GO TO 20
10 CONTINUE
20 C = DCOS(EA)
S = DSIN(EA)
B = DSQRT(1.0D0 - E * E)
RR = A * (1.0D0 - E * C)
VF = DSQRT(GME * A) / RR
DO 30 I = 1, 3
R(I) = A * ((C - E) * P(I) + B * S * Q(I))
V(I) = VF * (-S * P(I) + B * C * Q(I))
30 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C LMDESC: LM on the P64 approach. PMF (Moon centred, MF, km) is
C the commander's eye, the camera; body axes XB (thrust, up), YB
C (right), ZB (forward); LMALT the footpads' altitude (km).
C Range and footpad altitude fall as the square of the time to go
C to touchdown (the scenario's TOUCH event, LUT0), reaching 0 there;
C after it the LM stands landed. The line of sight to the site
C stays near 16 deg below the horizontal; the LM pitches up from
C 40 to 5 deg off vertical. This profile is ours, fitted to the
C film's descent frames (VIEW's pre-flight output). It is NOT the
C flown one: the altitude calls (Apollo Lunar Surface Journal,
C apollojournals.org/alsj/a11/a11.landing.html) have 1000 ft at
C 102:42:37, 300 at 102:43:46, 100 at 102:44:45, 40 at 102:45:17,
C 20 at 102:45:25 and "Contact Light" at 102:45:40; ours has
C about 3850, 1490, 350, 60 and 26 ft at those times. In the last
C minute ours falls about 350 ft where the flown descent fell
C about 100 ft, mostly a hover.
C Eye height above the footpads, LMEYE: ours, 5.1 m. The LM
C "stands 22 feet 11 inches high" with the gear out, the ascent
C stage is 12 feet 4 inches and the descent stage 10 feet 7 inches
C high (Apollo 11 press kit, printed pp. 96, 101), so the ascent
C stage's base is 3.23 m above the pads; we add 0.3 m to the cabin
C floor and 1.6 m for a standing man's eye (our estimates).
C-----------------------------------------------------------------------
SUBROUTINE LMDESC(GET, PMF, XB, YB, ZB)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, PMF(3), XB(3), YB(3), ZB(3)
DOUBLE PRECISION S(3), N(3), E(3), F(3), U(3), G(3)
DOUBLE PRECISION FI, LA, AZ, TAU, Q, X, H, TP, C, SN
INTEGER I
C The landing site and descent azimuth from the scenario.
FI = SNSLA(ISN) * DR
LA = SNSLO(ISN) * DR
AZ = SNSAZ(ISN) * DR
S(1) = DCOS(FI) * DCOS(LA)
S(2) = DCOS(FI) * DSIN(LA)
S(3) = DSIN(FI)
N(1) = -DSIN(FI) * DCOS(LA)
N(2) = -DSIN(FI) * DSIN(LA)
N(3) = DCOS(FI)
E(1) = -DSIN(LA)
E(2) = DCOS(LA)
E(3) = 0.0D0
DO 10 I = 1, 3
F(I) = DCOS(AZ) * N(I) + DSIN(AZ) * E(I)
10 CONTINUE
LMEYE = 5.1D-3
TAU = LUT0 - GET
IF (TAU .LT. 0.0D0) TAU = 0.0D0
IF (TAU .GT. 600.0D0) TAU = 600.0D0
Q = TAU / 250.0D0
X = 7.80D0 * Q * Q
H = 2.19D0 * Q * Q
LMALT = H
DO 20 I = 1, 3
G(I) = S(I) - (X / RM) * F(I)
20 CONTINUE
CALL VUNIT(G)
DO 30 I = 1, 3
U(I) = G(I)
30 CONTINUE
C Forward direction made horizontal at the LM.
C = F(1) * U(1) + F(2) * U(2) + F(3) * U(3)
DO 40 I = 1, 3
F(I) = F(I) - C * U(I)
40 CONTINUE
CALL VUNIT(F)
TP = (5.0D0 + 35.0D0 * Q) * DR
IF (TP .GT. 60.0D0 * DR) TP = 60.0D0 * DR
C = DCOS(TP)
SN = DSIN(TP)
DO 50 I = 1, 3
XB(I) = C * U(I) - SN * F(I)
ZB(I) = SN * U(I) + C * F(I)
C The eye, LMEYE above the footpads along the LM's up axis.
PMF(I) = (RM + H) * G(I) + LMEYE * XB(I)
50 CONTINUE
CALL VCRS(ZB, XB, YB)
RETURN
END
C
C-----------------------------------------------------------------------
C ERFIND: time TERISE at which the Earth's disc clears the lunar
C horizon, on the revolution ending at touchdown (the scenario's
C TOUCH event), or at its PHOTO event if it has no landing.
C-----------------------------------------------------------------------
SUBROUTINE ERFIND
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T, T1, T2, TM, F1, FM, P, ERCLR, TR, EVGET
INTEGER I, LEGAT
C The revolution ends at touchdown, or, in a scenario with no
C landing (Apollo 8), at its PHOTO event.
TR = LUT0
IF (TR .LE. 0.0D0) TR = EVGET(KEPHO)
P = 2.0D0 * PI / LGEL(11, LEGAT(TR, 2))
T1 = TR - P
F1 = ERCLR(T1)
DO 10 I = 1, 720
T = TR - P + DBLE(I) * P / 720.0D0
FM = ERCLR(T)
IF (F1 .LT. 0.0D0 .AND. FM .GE. 0.0D0) GO TO 20
T1 = T
F1 = FM
10 CONTINUE
TERISE = TR - 1800.0D0
RETURN
20 T2 = T
DO 30 I = 1, 40
TM = 0.5D0 * (T1 + T2)
FM = ERCLR(TM)
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (FM .LT. 0.0D0) THEN
T1 = TM
ELSE
T2 = TM
END IF
C RESTOMOD END
30 CONTINUE
TERISE = T2
RETURN
END
C
C ERCLR: angle (rad) of the Earth's lower limb above the lunar
C horizon as seen from the CSM.
DOUBLE PRECISION FUNCTION ERCLR(T)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION T, R(3), V(3), PM(3), E(3), DN(3), DE, VDOT
INTEGER I
CALL VSTATE(T, 2, R, V)
CALL MOONG(T, PM)
DO 10 I = 1, 3
E(I) = -PM(I) - R(I)
DN(I) = -R(I)
10 CONTINUE
DE = DSQRT(VDOT(E, E))
CALL VUNIT(E)
CALL VUNIT(DN)
ERCLR = DACOS(VDOT(E, DN)) - DASIN(RM / DSQRT(VDOT(R, R)))
& - DASIN(RE / DE)
RETURN
END
src/vlayer.f
C=======================================================================
C
C V I E W - 1 1 0 8 LAYER DISPATCHER
C
C Walks the scene's layer list and calls each layer by its id.
C Every layer is a subroutine of one file with the same argument
C list, (GET, VB, NV, SB, NS, LB, NL), and draws into the plot
C buffers in list order. FORTRAN 66 has no procedure variables,
C so the call is chosen by a computed GO TO over the id. One
C relocatable element of the kernel; see vdrive.f for the list.
C
C To add a layer: a new file with its subroutine, the next id, one
C GO TO target and CALL below, and the id in the scene lists.
C
C id layer element
C 1 plot frame and ticks lframe.f DFRAME
C 2 stars lstars.f DSTARS
C 3 Sun lsun.f DSUN
C 4 Moon and craters lmoon.f DMOON
C 5 Earth learth.f DEARTH
C 6 vehicles (placed models) lvehic.f MDRALL
C their labels and markers lvlab.f VLABEL
C 7 COAS reticle lcoas.f S7COAS
C 8 LM shadow lshad.f LMSHAD
C 9 LPD and LM window llpd.f OVLPD
C
C=======================================================================
SUBROUTINE LAYERS(GET, VB, NV, SB, NS, LB, NL)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, VB(5,MAXV), SB(3,MAXS), LB(4,MAXL)
INTEGER NV, NS, NL
INTEGER LL(12,9), K, L
C Layer lists, one column per scene, ended by 0. Scene 8 (the
C stack in translunar coast) has the sky and the vehicles. Scene 4
C (the LM pirouette) has no stars, as on the film (t22.png,
C t25.png). Scene 9 (Apollo 8 Earthrise) has scene 1's.
DATA LL / 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0,
& 1, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 8, 9, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0,
& 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 0 /
DO 90 K = 1, 12
L = LL(K, ISCN)
IF (L .LT. 1 .OR. L .GT. 9) GO TO 95
C Window overlays (COAS, LPD) only in the scene's window view.
IF (IVUSE .NE. 0 .AND. (L .EQ. 7 .OR. L .EQ. 9)) GO TO 90
GO TO (11, 12, 13, 14, 15, 16, 17, 18, 19), L
11 CALL DFRAME(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
12 CALL DSTARS(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
13 CALL DSUN(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
14 CALL DMOON(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
15 CALL DEARTH(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
16 CALL MDRALL(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
17 CALL S7COAS(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
18 CALL LMSHAD(GET, VB, NV, SB, NS, LB, NL)
GO TO 90
19 CALL OVLPD(GET, VB, NV, SB, NS, LB, NL)
90 CONTINUE
C The LM station view (in_view 3) carries the LM window overlay.
95 IF (IVUSE .EQ. 3) CALL OVLPD(GET, VB, NV, SB, NS, LB, NL)
RETURN
END
src/vmath.f
C=======================================================================
C
C V I E W - 1 1 0 8 VECTOR AND MATRIX
C
C Core element. 3-vector and 3x3 utilities. One relocatable
C element of the kernel; see vdrive.f for the list.
C
C=======================================================================
C
SUBROUTINE SETV(V, A, B, C)
DOUBLE PRECISION V(3), A, B, C
V(1) = A
V(2) = B
V(3) = C
RETURN
END
C
C=======================================================================
C VECTOR AND MATRIX UTILITIES
C=======================================================================
DOUBLE PRECISION FUNCTION VDOT(A, B)
DOUBLE PRECISION A(3), B(3)
VDOT = A(1) * B(1) + A(2) * B(2) + A(3) * B(3)
RETURN
END
C
DOUBLE PRECISION FUNCTION VNRM(A)
DOUBLE PRECISION A(3)
VNRM = DSQRT(A(1) * A(1) + A(2) * A(2) + A(3) * A(3))
RETURN
END
C
SUBROUTINE VCRS(A, B, C)
DOUBLE PRECISION A(3), B(3), C(3)
C(1) = A(2) * B(3) - A(3) * B(2)
C(2) = A(3) * B(1) - A(1) * B(3)
C(3) = A(1) * B(2) - A(2) * B(1)
RETURN
END
C
SUBROUTINE VUNIT(A)
DOUBLE PRECISION A(3), S
S = DSQRT(A(1) * A(1) + A(2) * A(2) + A(3) * A(3))
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (S .GT. 0.0D0) THEN
A(1) = A(1) / S
A(2) = A(2) / S
A(3) = A(3) / S
END IF
C RESTOMOD END
RETURN
END
C
C MXV: B = M A. MTXV: B = transpose(M) A. MXM: C = A B.
SUBROUTINE MXV(M, A, B)
DOUBLE PRECISION M(3,3), A(3), B(3)
INTEGER I
DO 10 I = 1, 3
B(I) = M(I,1) * A(1) + M(I,2) * A(2) + M(I,3) * A(3)
10 CONTINUE
RETURN
END
C
SUBROUTINE MTXV(M, A, B)
DOUBLE PRECISION M(3,3), A(3), B(3)
INTEGER I
DO 10 I = 1, 3
B(I) = M(1,I) * A(1) + M(2,I) * A(2) + M(3,I) * A(3)
10 CONTINUE
RETURN
END
C
SUBROUTINE MXM(A, B, C)
DOUBLE PRECISION A(3,3), B(3,3), C(3,3)
INTEGER I, J
DO 20 J = 1, 3
DO 10 I = 1, 3
C(I,J) = A(I,1) * B(1,J) + A(I,2) * B(2,J) + A(I,3) * B(3,J)
10 CONTINUE
20 CONTINUE
RETURN
END
C
C Active rotations by angle A (rad) about X, Y, Z.
SUBROUTINE ROTX(A, M)
DOUBLE PRECISION A, M(3,3)
CALL MUNIT(M)
M(2,2) = DCOS(A)
M(2,3) = -DSIN(A)
M(3,2) = DSIN(A)
M(3,3) = DCOS(A)
RETURN
END
C
SUBROUTINE ROTY(A, M)
DOUBLE PRECISION A, M(3,3)
CALL MUNIT(M)
M(1,1) = DCOS(A)
M(1,3) = DSIN(A)
M(3,1) = -DSIN(A)
M(3,3) = DCOS(A)
RETURN
END
C
SUBROUTINE ROTZ(A, M)
DOUBLE PRECISION A, M(3,3)
CALL MUNIT(M)
M(1,1) = DCOS(A)
M(1,2) = -DSIN(A)
M(2,1) = DSIN(A)
M(2,2) = DCOS(A)
RETURN
END
C
SUBROUTINE MUNIT(M)
DOUBLE PRECISION M(3,3)
INTEGER I, J
DO 20 J = 1, 3
DO 10 I = 1, 3
M(I,J) = 0.0D0
10 CONTINUE
M(J,J) = 1.0D0
20 CONTINUE
RETURN
END
src/vsrc.f
C=======================================================================
C
C V I E W - 1 1 0 8 STATE SOURCE
C
C The one place the display side gets the CSM's state. It reads
C the replay (the scenario's legs, traj.f) or, when the frame asks
C for it (in_src, which the chassis passes as in_flags bit 3) and
C the tape covers the time, the tape
C the engine wrote (tape.f, sim.f). Scene cameras and layers do
C not know which. One relocatable element of the kernel; see
C vdrive.f for the list.
C
C=======================================================================
C
C VSTATE: the CSM at GET, about the Earth (IBODY = 1, geocentric
C EQ) or the Moon (IBODY = 2, selenocentric EQ), km and km/s.
C ISRCU records the source used: 0 replay, 1 sim with state
C vector updates, 2 sim without.
SUBROUTINE VSTATE(GET, IBODY, R, V)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, R(3), V(3), PM(3), VM(3)
INTEGER IBODY, IOK, I
ISRCU = 0
IF (ISRC .NE. 1 .OR. ITPSN .NE. ISN) GO TO 20
CALL TPGET(1, GET, R, V, IOK)
IF (IOK .EQ. 0) GO TO 20
ISRCU = 2 - MOD(ISIMF, 2)
IF (IBODY .EQ. 1) RETURN
CALL MOONV(GET, PM, VM)
DO 10 I = 1, 3
R(I) = R(I) - PM(I)
V(I) = V(I) - VM(I)
10 CONTINUE
RETURN
20 IF (IBODY .EQ. 1) CALL ERTORB(GET, R, V)
IF (IBODY .EQ. 2) CALL LUNORB(GET, R, V)
RETURN
END
C
C SIMERR: for the header, the reference row nearest GET that the
C last engine run passed: its g.e.t. and the position (km) and
C velocity (ft/s) error there. JN = 0 if there is none.
SUBROUTINE SIMERR(GET, JN, TR, EP, EV)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, TR, EP, EV, D, DB
INTEGER JN, J
JN = 0
TR = 0.0D0
EP = 0.0D0
EV = 0.0D0
IF (ITPSN .NE. ISN) RETURN
DB = 1.0D30
DO 10 J = 1, NREF
IF (RFOK(J) .EQ. 0) GO TO 10
D = DABS(RFP(3,J) - GET)
IF (D .GE. DB) GO TO 10
DB = D
JN = J
10 CONTINUE
IF (JN .EQ. 0) RETURN
TR = RFP(3,JN)
EP = RFERR(1,JN)
EV = RFERR(2,JN)
RETURN
END
src/vtext.f
C=======================================================================
C
C V I E W - 1 1 0 8 TEXT
C
C Core element. Records for the recorder's character
C generator. One relocatable element of
C the kernel; see vdrive.f for the list.
C
C=======================================================================
C
C=======================================================================
C TEXT. Records for the recorder's character generator (the SC-4020
C class had a "type character" order; docs/univac-1108.md): X, Y
C of the first character's lower left (plot deg), height (plot
C deg), start index in TC; each string is character codes ending
C in 0. Tick numbers when IFLG bit 1 is set, at the ticks DFRAME
C draws, OUTSIDE the box: left edge right-aligned, right edge,
C bottom edge centred below, 1.0 percent of the field high (read
C from the film, descent_t29.png, and the report's plot pages).
C Names of nav stars, Sun, Earth and Moon beside their labels when
C bit 0 is set, 1.4 percent of the field high, by the label level
C (ILEV): primary Sun, Earth, Moon, the landing site, vehicles and
C the pad; secondary also nav stars and maria; all, everything.
C The labels (LB) themselves are not filtered. Also vehicles (LB
C kind 8: CM, SM, LM, S-IVB, CSM by id) and the launch pad (kind 9,
C its name from the PAD card), both modern additions. Crater
C names stay with the page (LB kind 2). Character width 0.7 of
C the height, used for alignment: ours.
C=======================================================================
SUBROUTINE TXALL(LB, NL, TB, NT, TC, NCH)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION LB(4,MAXL), TB(4,MAXT)
INTEGER NL, NT, TC(MAXTC), NCH
DOUBLE PRECISION B, ST, TL, H, V
INTEGER I, J, K, N, ID, IC(24), VEHCH(25)
C Vehicle names, 5 codes each, zero padded: CM, SM, LM, S-IVB, CSM.
DATA VEHCH / 67, 77, 0, 0, 0, 83, 77, 0, 0, 0, 76, 77, 0, 0, 0,
& 83, 45, 73, 86, 66, 67, 83, 77, 0, 0 /
NT = 0
NCH = 0
B = BOXH
H = 0.02D0 * B
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (MOD(IFLG / 2, 2) .EQ. 1) THEN
ST = 20.0D0
IF (2.0D0 * B .LE. 60.0D0) ST = 10.0D0
IF (2.0D0 * B .LE. 25.0D0) ST = 5.0D0
TL = 0.02D0 * B
N = INT(B / ST + 1.0D-9)
DO 10 K = -N, N
V = DBLE(K) * ST
IF (DABS(V) .GE. B - 1.0D-9) GO TO 10
CALL ITOC(NINT(V), IC, J)
CALL TXPUT(TB, NT, TC, NCH, V - 0.35D0 * H * DBLE(J),
& -B - 1.5D0 * H, H, IC, J)
CALL TXPUT(TB, NT, TC, NCH, -B - 0.5D0 * H
& - 0.7D0 * H * DBLE(J), V - 0.5D0 * H, H, IC, J)
CALL TXPUT(TB, NT, TC, NCH, B + 0.5D0 * H,
& V - 0.5D0 * H, H, IC, J)
10 CONTINUE
END IF
C RESTOMOD END
IF (MOD(IFLG, 2) .EQ. 0) RETURN
H = 0.028D0 * B
DO 40 I = 1, NL
K = NINT(LB(3,I))
ID = NINT(LB(4,I))
N = 0
IF ((K .EQ. 1 .OR. K .EQ. 6) .AND. ILEV .LT. 2) GO TO 40
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (K .EQ. 1 .AND. ID .GE. 1 .AND. ID .LE. NNAV) THEN
DO 20 J = 1, 10
IF (NAVCH((ID - 1) * 10 + J) .EQ. 0) GO TO 30
N = N + 1
IC(N) = NAVCH((ID - 1) * 10 + J)
20 CONTINUE
ELSE IF (K .GE. 3 .AND. K .LE. 5) THEN
DO 25 J = 1, 5
IF (BODCH((K - 3) * 5 + J) .EQ. 0) GO TO 30
N = N + 1
IC(N) = BODCH((K - 3) * 5 + J)
25 CONTINUE
ELSE IF (K .EQ. 6 .AND. ID .GE. 1 .AND. ID .LE. NMARE) THEN
C Mare names centred on the mare's centre.
DO 26 J = 1, 24
IF (MRCH((ID - 1) * 24 + J) .EQ. 0) GO TO 27
N = N + 1
IC(N) = MRCH((ID - 1) * 24 + J)
26 CONTINUE
27 IF (N .GT. 0) CALL TXPUT(TB, NT, TC, NCH,
& LB(1,I) - 0.35D0 * H * DBLE(N), LB(2,I) - 0.5D0 * H,
& H, IC, N)
GO TO 40
ELSE IF (K .EQ. 7) THEN
DO 28 J = 1, 22
N = N + 1
IC(N) = SITECH(J)
28 CONTINUE
ELSE IF (K .EQ. 8 .AND. ID .GE. 1 .AND. ID .LE. 5) THEN
DO 29 J = 1, 5
IF (VEHCH((ID - 1) * 5 + J) .EQ. 0) GO TO 30
N = N + 1
IC(N) = VEHCH((ID - 1) * 5 + J)
29 CONTINUE
ELSE IF (K .EQ. 9) THEN
DO 32 J = 1, 8
IF (PADCH(8 * (ISN - 1) + J) .EQ. 0) GO TO 30
N = N + 1
IC(N) = PADCH(8 * (ISN - 1) + J)
32 CONTINUE
END IF
C RESTOMOD END
30 IF (N .GT. 0) CALL TXPUT(TB, NT, TC, NCH, LB(1,I) + 0.4D0 * H,
& LB(2,I) + 0.4D0 * H, H, IC, N)
40 CONTINUE
RETURN
END
C
C TXPUT: append a text record of N codes IC.
SUBROUTINE TXPUT(TB, NT, TC, NCH, X, Y, H, IC, N)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION TB(4,MAXT), X, Y, H
INTEGER NT, TC(MAXTC), NCH, IC(24), N, K
IF (NT .GE. MAXT .OR. NCH + N + 1 .GT. MAXTC) RETURN
NT = NT + 1
TB(1,NT) = X
TB(2,NT) = Y
TB(3,NT) = H
TB(4,NT) = DBLE(NCH + 1)
DO 10 K = 1, N
TC(NCH + K) = IC(K)
10 CONTINUE
NCH = NCH + N + 1
TC(NCH) = 0
RETURN
END
C
C ITOC: integer IV to character codes IC(1..N), ASCII digits.
SUBROUTINE ITOC(IV, IC, N)
INTEGER IV, IC(24), N, M, K, D(10), ND
M = IABS(IV)
ND = 0
10 ND = ND + 1
D(ND) = MOD(M, 10)
M = M / 10
IF (M .GT. 0 .AND. ND .LT. 10) GO TO 10
N = 0
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IV .LT. 0) THEN
N = 1
IC(1) = 45
END IF
C RESTOMOD END
DO 20 K = ND, 1, -1
N = N + 1
IC(N) = 48 + D(K)
20 CONTINUE
RETURN
END
src/vview.f
C=======================================================================
C
C V I E W - 1 1 0 8 CAMERA POINTING
C
C Core element. The camera target and the external view, applied
C after the scene has set its own camera (SCNCAM) and placed its
C models (SCNMOD). One relocatable element of the kernel; see
C vdrive.f for the list.
C
C TN D-6853 (printed p. 13) offers "an inertially fixed platform
C or a local-vertical platform". Pointing at a target and flying
C the camera around one are a third and fourth mode: ours.
C Window view (in_view 0) with a target: the reference boresight
C points from the scene's camera at the target; free-look
C yaw, pitch and roll are offsets from it.
C Stations (in_view 2, 3): the camera at the CM or LM eye, the
C cabin around it (STATCM, STATLM).
C External view (in_view 1): the camera sits D from the target
C and looks at it; yaw and pitch carry it around the target
C (starting from the side the scene's own camera is on), roll
C turns the picture, the field of view zooms. D is ours: 60 m
C for the CSM or the docked stack, 40 m for the LM, 60,000 km
C for the Earth, 36,737 km (35,000 km up, as scene 6) for the
C Moon. The Sun is a target only for the window view.
C Scene 6 keeps its own Moon-centred camera and ignores both.
C=======================================================================
C
C-----------------------------------------------------------------------
C VIEWPT: apply the view and target to the camera: CG (geocentric,
C km) may move, BREF, UREF, RREF may turn, and placed models are
C carried along. LOOKD = 1 if the camera axes are already set
C (external view), so VFRAME skips its free-look step.
C-----------------------------------------------------------------------
SUBROUTINE VIEWPT(GET, PM, CG, YAW, PIT, ROL, LOOKD)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, PM(3), CG(3), YAW, PIT, ROL
INTEGER LOOKD
DOUBLE PRECISION TG(3), D(3), DIST, DS(3), CG0(3), P(3), VDOT
DOUBLE PRECISION C
INTEGER IT, IOK, I
LOOKD = 0
IVUSE = IVIEW
C Scene 8's own view is the external one.
IF (ISCN .EQ. 8 .AND. IVUSE .EQ. 0) IVUSE = 1
C Scene 5 is the LM station already.
IF (ISCN .EQ. 5 .AND. IVUSE .EQ. 3) IVUSE = 0
IF (ISCN .EQ. 6) IVUSE = 0
IT = ITARG
IF (ISCN .EQ. 6) RETURN
C Stations: the camera moves to the eye, the cabin goes around it.
C Where the scene has no such vehicle, the window view.
IF (IVUSE .EQ. 2) CALL STATCM(CG, IOK)
IF (IVUSE .EQ. 2 .AND. IOK .EQ. 0) IVUSE = 0
IF (IVUSE .EQ. 3) CALL STATLM(CG, IOK)
IF (IVUSE .EQ. 3 .AND. IOK .EQ. 0) IVUSE = 0
C The LM station keeps its window's own aim (its overlay is drawn in
C the reference frame, OVLPD).
IF (IVUSE .EQ. 3) RETURN
IF (IVUSE .NE. 1 .AND. IT .EQ. 0) RETURN
C The external view's default target is the scene's subject.
IF (IVUSE .EQ. 1 .AND. (IT .EQ. 0 .OR. IT .EQ. 3))
& CALL TGTDEF(IT)
C Seen from outside, a camera riding the CSM (scenes 1, 2, 3, 4,
C 7, 9) shows the CSM: its outline with the CM's base 1.2 m behind
C the eye along the scene's boresight, its X axis along that
C boresight (ours).
IF (IVUSE .EQ. 1 .AND. MDON(KCSM) .EQ. 0 .AND. ISCN .NE. 5
& .AND. ISCN .NE. 8) CALL CSMCAM
CALL TGTPOS(GET, IT, PM, CG, TG, DIST, IOK)
IF (IOK .EQ. 0) RETURN
DO 10 I = 1, 3
CG0(I) = CG(I)
D(I) = TG(I) - CG(I)
10 CONTINUE
C Reference boresight at the target, up kept as near the scene's
C as it can be.
IF (VDOT(D, D) .LE. 1.0D-18) GO TO 22
CALL VUNIT(D)
DO 20 I = 1, 3
BREF(I) = D(I)
20 CONTINUE
22 CONTINUE
C = VDOT(UREF, BREF)
DO 25 I = 1, 3
UREF(I) = UREF(I) - C * BREF(I)
25 CONTINUE
IF (VDOT(UREF, UREF) .LT. 1.0D-12) CALL PERP(BREF, UREF, P)
CALL VUNIT(UREF)
CALL VCRS(BREF, UREF, RREF)
CALL VUNIT(RREF)
IF (IVUSE .NE. 1) RETURN
C
C External: free-look sets the direction, the camera backs off
C along it to DIST from the target.
CALL LOOK(YAW, PIT, ROL)
LOOKD = 1
DO 30 I = 1, 3
CG(I) = TG(I) - DIST * CB(I)
DS(I) = CG0(I) - CG(I)
30 CONTINUE
C Carry the placed models, shifted by DS.
CALL MSHIFT(DS, 0)
RETURN
END
C
C TGTDEF: the scene's default external target: the reference body
C (Moon in scenes 1, 5, 6, 9; Earth in 2, 3), the LM (4, 7), the
C CSM (8).
SUBROUTINE TGTDEF(IT)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
INTEGER IT, ID(9)
DATA ID / 2, 1, 1, 5, 2, 2, 5, 4, 2 /
IT = ID(ISCN)
RETURN
END
C
C-----------------------------------------------------------------------
C TGTPOS: target IT's geocentric position TG (km) and the external
C view's distance DIST (km). IOK = 0 if the scene has no such
C target, or it is the camera itself in a window view. Vehicles:
C the placed models (the CSM, KCSM; the LM, KLMD or KLMS), else
C the CSM from the state source (VSTATE) where the camera is not
C the CSM (scenes 5, 6).
C-----------------------------------------------------------------------
SUBROUTINE TGTPOS(GET, IT, PM, CG, TG, DIST, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION GET, PM(3), CG(3), TG(3), DIST, R(3), V(3)
DOUBLE PRECISION D
INTEGER IT, IOK, I, KL
IOK = 1
DIST = 0.0D0
DO 5 I = 1, 3
TG(I) = 0.0D0
5 CONTINUE
C RESTOMOD BEGIN: block IF is FORTRAN 77 (1978)
IF (IT .EQ. 1) THEN
DIST = 60000.0D0
ELSE IF (IT .EQ. 2) THEN
DO 10 I = 1, 3
TG(I) = PM(I)
10 CONTINUE
DIST = RM + 35000.0D0
ELSE IF (IT .EQ. 3) THEN
DO 20 I = 1, 3
TG(I) = CG(I) + 1.495978707D8 * SUNU(I)
20 CONTINUE
IF (IVUSE .EQ. 1) IOK = 0
ELSE IF (IT .EQ. 4) THEN
DIST = 0.060D0
IF (MDON(KCSM) .EQ. 1) THEN
C The placed CSM: aim at the top of its tunnel.
DO 30 I = 1, 3
TG(I) = CG(I) + MDP(I,KCSM) + 3.2D-3 * MDAT(I,1,KCSM)
30 CONTINUE
ELSE IF (ISCN .EQ. 5 .OR. ISCN .EQ. 6) THEN
CALL VSTATE(GET, 1, R, V)
DO 35 I = 1, 3
TG(I) = R(I)
35 CONTINUE
ELSE
IOK = 0
END IF
ELSE
DIST = 0.040D0
KL = 0
IF (MDON(KLMD) .EQ. 1) KL = KLMD
IF (MDON(KLMS) .EQ. 1) KL = KLMS
IF (KL .EQ. 0) THEN
IOK = 0
ELSE
DO 50 I = 1, 3
TG(I) = CG(I) + MDP(I,KL)
50 CONTINUE
END IF
END IF
C RESTOMOD END
IF (IOK .EQ. 0 .OR. IVUSE .EQ. 1) RETURN
C A window view cannot aim at its own camera.
D = 0.0D0
DO 60 I = 1, 3
D = D + (TG(I) - CG(I))**2
60 CONTINUE
IF (D .LT. 1.0D-10) IOK = 0
RETURN
END
C
C CSMCAM: place the CSM outline around the scene's camera, the CM's
C base 1.2 m behind the eye, X along the reference boresight, Z
C along its up (ours).
SUBROUTINE CSMCAM
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION AT(3,3), P(3), Z(3)
INTEGER I
DO 10 I = 1, 3
AT(I,1) = BREF(I)
AT(I,3) = UREF(I)
P(I) = -1.2D-3 * BREF(I)
Z(I) = 0.0D0
10 CONTINUE
CALL VCRS(AT(1,3), AT(1,1), AT(1,2))
CALL MPLACE(KCSM, AT, P, Z)
RETURN
END
C
C MSHIFT: place every placed model again, shifted by DS (km); model
C KX (0 for none) is left out, as the one the camera is inside.
SUBROUTINE MSHIFT(DS, KX)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION DS(3), AT(3,3), P(3), BO(3)
INTEGER KX, KP(MMOD), I, J, K
DO 10 K = 1, NMOD
KP(K) = MDON(K)
10 CONTINUE
CALL MCLEAR
DO 30 K = 1, NMOD
IF (KP(K) .EQ. 0 .OR. K .EQ. KX) GO TO 30
DO 20 I = 1, 3
P(I) = MDP(I,K) + DS(I)
BO(I) = MDBO(I,K)
DO 15 J = 1, 3
AT(I,J) = MDAT(I,J,K)
15 CONTINUE
20 CONTINUE
CALL MPLACE(K, AT, P, BO)
30 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C STATCM: the CM station. The camera at the CM eye (CMEYE), looking
C along the CSM's +X axis with its -Z up, the view of MSC IN 69-FM-
C 197's CSM maneuver plots (see CMCAB); the cabin (KCMC) around it.
C The CSM's axes: the placed CSM's (scene 8), else X along the
C scene's boresight and Z against its up, so the station looks
C where the window view looked (ours). Not in scenes 5 and 6.
C-----------------------------------------------------------------------
SUBROUTINE STATCM(CG, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION CG(3), AT(3,3), E(3), V(3), W(3), DS(3), Z(3)
INTEGER IOK, I, J
IOK = 0
IF (ISCN .EQ. 5 .OR. ISCN .EQ. 6) RETURN
IOK = 1
CALL CMEYE(E)
CALL SETV(Z, 0.0D0, 0.0D0, 0.0D0)
IF (MDON(KCSM) .EQ. 0) GO TO 20
DO 10 I = 1, 3
V(I) = (E(I) - MDBO(I,KCSM)) * 1.0D-3
DO 5 J = 1, 3
AT(I,J) = MDAT(I,J,KCSM)
5 CONTINUE
10 CONTINUE
CALL MXV(AT, V, W)
DO 15 I = 1, 3
W(I) = W(I) + MDP(I,KCSM)
CG(I) = CG(I) + W(I)
DS(I) = -W(I)
15 CONTINUE
CALL MSHIFT(DS, KCSM)
GO TO 30
20 DO 25 I = 1, 3
AT(I,1) = BREF(I)
AT(I,3) = -UREF(I)
25 CONTINUE
CALL VCRS(AT(1,3), AT(1,1), AT(1,2))
30 DO 35 I = 1, 3
BREF(I) = AT(I,1)
UREF(I) = -AT(I,3)
35 CONTINUE
CALL VCRS(BREF, UREF, RREF)
CALL MPLACE(KCMC, AT, Z, E)
RETURN
END
C
C-----------------------------------------------------------------------
C STATLM: the LM station. The camera at the commander's eye in the
C placed LM (scenes 4, 7, 8), 3.2 m up its X axis, 0.5 m to -Y and
C 0.8 m forward (ours), looking as scene 5 does, 46 deg down from
C the LM's +Z in the X-Z plane; the LM window and LPD overlay
C (OVLPD) is drawn about that aim. Not where no LM is placed.
C-----------------------------------------------------------------------
SUBROUTINE STATLM(CG, IOK)
C RESTOMOD BEGIN: file INCLUDE; FORTRAN V's named PDP elements
INCLUDE 'viewdims.inc'
INCLUDE 'viewcom.inc'
C RESTOMOD END
DOUBLE PRECISION CG(3), AT(3,3), E(3), V(3), W(3), DS(3), C, S
INTEGER IOK, I, J, KL
IOK = 0
KL = 0
IF (MDON(KLMD) .EQ. 1) KL = KLMD
IF (MDON(KLMS) .EQ. 1) KL = KLMS
IF (KL .EQ. 0) RETURN
IOK = 1
CALL SETV(E, 3.2D0, -0.5D0, 0.8D0)
DO 10 I = 1, 3
V(I) = (E(I) - MDBO(I,KL)) * 1.0D-3
DO 5 J = 1, 3
AT(I,J) = MDAT(I,J,KL)
5 CONTINUE
10 CONTINUE
CALL MXV(AT, V, W)
DO 15 I = 1, 3
W(I) = W(I) + MDP(I,KL)
CG(I) = CG(I) + W(I)
DS(I) = -W(I)
15 CONTINUE
CALL MSHIFT(DS, KL)
C = DCOS(46.0D0 * DR)
S = DSIN(46.0D0 * DR)
DO 20 I = 1, 3
BREF(I) = C * AT(I,3) - S * AT(I,1)
UREF(I) = S * AT(I,3) + C * AT(I,1)
20 CONTINUE
CALL VCRS(BREF, UREF, RREF)
CALL VUNIT(RREF)
RETURN
END