1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
//
// GENERATED FILE
//
use super::*;
use f2rust_std::*;
pub const UBPL: i32 = 4;
struct SaveVars {
Z: StackArray<f64, 3>,
}
impl SaveInit for SaveVars {
fn new() -> Self {
let mut Z = StackArray::<f64, 3>::new(1..=3);
{
use f2rust_std::data::Val;
let mut clist = [Val::D(0.0), Val::D(0.0), Val::D(1.0)].into_iter();
Z.iter_mut()
.for_each(|n| *n = clist.next().unwrap().into_f64());
debug_assert!(clist.next().is_none(), "DATA not fully initialised");
}
Self { Z }
}
}
//$Procedure ZZSGLATX ( Line segment latitude extent )
pub fn ZZSGLATX(
P1: &[f64],
P2: &[f64],
MINLAT: &mut f64,
MINP: &mut [f64],
MAXLAT: &mut f64,
MAXP: &mut [f64],
ctx: &mut Context,
) -> f2rust_std::Result<()> {
let save = ctx.get_vars::<SaveVars>();
let save = &mut *save.borrow_mut();
let P1 = DummyArray::new(P1, 1..=3);
let P2 = DummyArray::new(P2, 1..=3);
let mut MINP = DummyArrayMut::new(MINP, 1..=3);
let mut MAXP = DummyArrayMut::new(MAXP, 1..=3);
let mut CREASE = StackArray::<f64, 3>::new(1..=3);
let mut DIR = StackArray::<f64, 3>::new(1..=3);
let mut DP1: f64 = 0.0;
let mut DP2: f64 = 0.0;
let mut LAT: f64 = 0.0;
let mut LAT1: f64 = 0.0;
let mut LAT2: f64 = 0.0;
let mut LON: f64 = 0.0;
let mut NORMAL = StackArray::<f64, 3>::new(1..=3);
let mut PLANE2 = StackArray::<f64, 4>::new(1..=UBPL);
let mut R: f64 = 0.0;
let mut T = StackArray::<f64, 3>::new(1..=3);
let mut NXPTS: i32 = 0;
//
// SPICELIB functions
//
//
// Local variables
//
//
// Saved variables
//
//
// Initial values
//
//
// Standard SPICE error handling.
//
if RETURN(ctx) {
return Ok(());
}
CHKIN(b"ZZSGLATX", ctx)?;
//
// Start by computing latitude at the segment's endpoints.
//
RECLAT(P1.as_slice(), &mut R, &mut LON, &mut LAT1);
RECLAT(P2.as_slice(), &mut R, &mut LON, &mut LAT2);
//
// Initialize the outputs using latitudes of the endpoints.
// If there are interior extrema, we'll update these outputs
// as needed.
//
if (LAT1 <= LAT2) {
*MINLAT = LAT1;
*MAXLAT = LAT2;
VEQU(P1.as_slice(), MINP.as_slice_mut());
VEQU(P2.as_slice(), MAXP.as_slice_mut());
} else {
*MINLAT = LAT2;
*MAXLAT = LAT1;
VEQU(P2.as_slice(), MINP.as_slice_mut());
VEQU(P1.as_slice(), MAXP.as_slice_mut());
}
//
// We want to work with the plane containing the origin, P1, and P2.
// We'll call this plane PLANE1. First see whether P1 and P2 are
// linearly independent.
//
VCRSS(P1.as_slice(), P2.as_slice(), NORMAL.as_slice_mut());
if VZERO(NORMAL.as_slice()) {
//
// We have a special case: P1 and P2 define a line passing
// through the origin. The latitude extrema lie on the
// segment endpoints, and possibly at every point on the
// segment. We've already computed the outputs.
//
CHKOUT(b"ZZSGLATX", ctx)?;
return Ok(());
}
//
// At this point we know that NORMAL is non-zero. Convert it
// to a unit vector.
//
VHATIP(NORMAL.as_slice_mut());
//
// Let ALPHA be the non-negative angle between PLANE1 and the X-Y
// plane. Then ALPHA and -ALPHA are, respectively, the maximum and
// minimum possible latitudes attained on the input segment.
// However, these values are not necessarily attained on the
// segment; we'll need to perform further analysis to find out. We
// don't need to compute ALPHA, but we'll refer to it in the
// discussion below.
//
// The next step is to find the normal vector to the plane defined
// by Z and NORMAL. We'll call this plane PLANE2. This plane might
// not exist if NORMAL and Z are linearly dependent. If PLANE2
// does exist, the X-Y plane and PLANE1 intersect in a "crease"
// that is normal to PLANE2.
//
VCRSS(save.Z.as_slice(), NORMAL.as_slice(), CREASE.as_slice_mut());
if VZERO(CREASE.as_slice()) {
//
// Z and NORMAL are linearly dependent; PLANE1 coincides (up to
// round-off error) with the X-Y plane. We've already computed
// the outputs.
//
CHKOUT(b"ZZSGLATX", ctx)?;
return Ok(());
}
//
// At this point we know CREASE is non-zero. Convert
// it to a unit vector.
//
VHATIP(CREASE.as_slice_mut());
//
// By construction, CREASE is orthogonal to NORMAL. PLANE2
// cuts PLANE1 in a line L passing through the origin. If
// the line segment has an interior latitude extremum,
// the point T where that extremum is attained lies on L.
// The segment is tangent at T to a nappe of a cone, centered on
// the Z-axis and having its apex at the origin, for which
// the half-angle is (pi/2)-ALPHA. The point T lies in PLANE2
// since L is contained in PLANE2.
//
// If a single tangent point T exists in the interior of the
// segment, then the endpoints must be on opposite sides of PLANE2.
// See whether this is the case.
//
DP1 = VDOT(P1.as_slice(), CREASE.as_slice());
DP2 = VDOT(P2.as_slice(), CREASE.as_slice());
if OPSGND(DP1, DP2) {
//
// The segment crosses PLANE2 at an interior point; this
// point is where the extremum occurs. Solve for the
// intersection.
//
// CREASE is guaranteed to be a unit vector. A zero input
// vector is the only cause for which NVC2PL will signal
// an error. Therefore we don't check FAILED after the
// following call.
//
NVC2PL(CREASE.as_slice(), 0.0, PLANE2.as_slice_mut(), ctx)?;
VSUB(P2.as_slice(), P1.as_slice(), DIR.as_slice_mut());
INRYPL(
P1.as_slice(),
DIR.as_slice(),
PLANE2.as_slice(),
&mut NXPTS,
T.as_slice_mut(),
ctx,
)?;
if FAILED(ctx) {
CHKOUT(b"ZZSGLATX", ctx)?;
return Ok(());
}
if (NXPTS == 1) {
//
// This is the normal case: we have one intersection of the
// segment with PLANE2. Update whichever of the extrema is
// superseded by the interior value.
//
// Note that this case can occur when NORMAL is orthogonal to
// Z, making ALPHA equal to pi/2. The nappes are degenerate in
// this case, consisting of the positive and negative Z-axes.
// This degenerate case occurs when the segment intersects the
// Z-axis in a point other than the origin, and the endpoints
// are linearly independent.
//
// This is not a special case computationally.
//
RECLAT(T.as_slice(), &mut R, &mut LON, &mut LAT);
if (LAT > *MAXLAT) {
*MAXLAT = LAT;
VEQU(T.as_slice(), MAXP.as_slice_mut());
} else if (LAT < *MINLAT) {
*MINLAT = LAT;
VEQU(T.as_slice(), MINP.as_slice_mut());
}
//
// There can be only one local extremum, so we're done.
//
}
//
// If NXPTS is not 1, then even though the endpoints are on
// opposite sides of PLANE2, either the segment was found to lie
// in PLANE2 or no intersection was found. This situation must be
// due to finite precision arithmetic error. We'll make do with
// the extrema already found.
}
//
// We reach this point if we found a local extremum or if any of the
// following are true:
//
// 1) The segment misses PLANE2 altogether, in which case
// there's no tangency point.
//
// 2) One endpoint lies on PLANE2 and one endpoint does not.
//
// 3) Both endpoints lie in PLANE2. Then both endpoints lie
// in L, so we should have found them to be linearly
// dependent. This situation must be due to finite precision
// arithmetic error.
//
// In all of the numbered cases the extrema occur at the endpoints.
// and have been found already. In all cases, the outputs are set.
//
CHKOUT(b"ZZSGLATX", ctx)?;
Ok(())
}