rsspice/generated/spicelib/dskw02.rs
1//
2// GENERATED FILE
3//
4
5use super::*;
6use crate::SpiceContext;
7use f2rust_std::*;
8
9const SRFIDX: i32 = 1;
10const CTRIDX: i32 = (SRFIDX + 1);
11const CLSIDX: i32 = (CTRIDX + 1);
12const TYPIDX: i32 = (CLSIDX + 1);
13const FRMIDX: i32 = (TYPIDX + 1);
14const SYSIDX: i32 = (FRMIDX + 1);
15const PARIDX: i32 = (SYSIDX + 1);
16const NSYPAR: i32 = 10;
17const MN1IDX: i32 = (PARIDX + NSYPAR);
18const MX1IDX: i32 = (MN1IDX + 1);
19const MN2IDX: i32 = (MX1IDX + 1);
20const MX2IDX: i32 = (MN2IDX + 1);
21const MN3IDX: i32 = (MX2IDX + 1);
22const MX3IDX: i32 = (MN3IDX + 1);
23const BTMIDX: i32 = (MX3IDX + 1);
24const ETMIDX: i32 = (BTMIDX + 1);
25const DSKDSZ: i32 = ETMIDX;
26const SVFCLS: i32 = 1;
27const GENCLS: i32 = 2;
28const LATSYS: i32 = 1;
29const CYLSYS: i32 = 2;
30const RECSYS: i32 = 3;
31const PDTSYS: i32 = 4;
32const XFRACT: f64 = 0.0000000001;
33const KEYXFR: i32 = 1;
34const SGREED: f64 = 0.00000001;
35const KEYSGR: i32 = (KEYXFR + 1);
36const SGPADM: f64 = 0.0000000001;
37const KEYSPM: i32 = (KEYSGR + 1);
38const PTMEMM: f64 = 0.0000001;
39const KEYPTM: i32 = (KEYSPM + 1);
40const ANGMRG: f64 = 0.000000000001;
41const KEYAMG: i32 = (KEYPTM + 1);
42const LONALI: f64 = 0.000000000001;
43const KEYLAL: i32 = (KEYAMG + 1);
44const IXNV: i32 = 1;
45const IXNP: i32 = (IXNV + 1);
46const IXNVXT: i32 = (IXNP + 1);
47const IXVGRX: i32 = (IXNVXT + 1);
48const IXCGSC: i32 = (IXVGRX + 3);
49const IXVXPS: i32 = (IXCGSC + 1);
50const IXVXLS: i32 = (IXVXPS + 1);
51const IXVTLS: i32 = (IXVXLS + 1);
52const IXPLAT: i32 = (IXVTLS + 1);
53const IXDSCR: i32 = 1;
54const DSCSZ2: i32 = 24;
55const IXVTBD: i32 = (IXDSCR + DSCSZ2);
56const IXVXOR: i32 = (IXVTBD + 6);
57const IXVXSZ: i32 = (IXVXOR + 3);
58const IXVERT: i32 = (IXVXSZ + 1);
59const KWNV: i32 = 1;
60const KWNP: i32 = (KWNV + 1);
61const KWNVXT: i32 = (KWNP + 1);
62const KWVGRX: i32 = (KWNVXT + 1);
63const KWCGSC: i32 = (KWVGRX + 1);
64const KWVXPS: i32 = (KWCGSC + 1);
65const KWVXLS: i32 = (KWVXPS + 1);
66const KWVTLS: i32 = (KWVXLS + 1);
67const KWPLAT: i32 = (KWVTLS + 1);
68const KWVXPT: i32 = (KWPLAT + 1);
69const KWVXPL: i32 = (KWVXPT + 1);
70const KWVTPT: i32 = (KWVXPL + 1);
71const KWVTPL: i32 = (KWVTPT + 1);
72const KWCGPT: i32 = (KWVTPL + 1);
73const KWDSC: i32 = (KWCGPT + 1);
74const KWVTBD: i32 = (KWDSC + 1);
75const KWVXOR: i32 = (KWVTBD + 1);
76const KWVXSZ: i32 = (KWVXOR + 1);
77const KWVERT: i32 = (KWVXSZ + 1);
78const MAXVRT: i32 = 16000002;
79const MAXPLT: i32 = (2 * (MAXVRT - 2));
80const MAXNPV: i32 = (((3 * MAXPLT) / 2) + 1);
81const MAXVOX: i32 = 100000000;
82const MAXCGR: i32 = 100000;
83const MAXEDG: i32 = 120;
84const SIVGRX: i32 = 1;
85const SICGSC: i32 = (SIVGRX + 3);
86const SIVXNP: i32 = (SICGSC + 1);
87const SIVXNL: i32 = (SIVXNP + 1);
88const SIVTNL: i32 = (SIVXNL + 1);
89const SICGRD: i32 = (SIVTNL + 1);
90const IXIFIX: i32 = (MAXCGR + 7);
91const SIVTBD: i32 = 1;
92const SIVXOR: i32 = (SIVTBD + 6);
93const SIVXSZ: i32 = (SIVXOR + 3);
94const IXDFIX: i32 = 10;
95const MAXVXP: i32 = (MAXPLT / 2);
96const MAXCEL: i32 = 60000000;
97const MXNVLS: i32 = (MAXCEL + (MAXVXP / 2));
98const SPAISZ: i32 = ((((IXIFIX + MAXVXP) + MXNVLS) + MAXVRT) + MAXNPV);
99
100/// DSK, write type 2 segment
101///
102/// Write a type 2 segment to a DSK file.
103///
104/// # Required Reading
105///
106/// * [DAS](crate::required_reading::das)
107/// * [DSK](crate::required_reading::dsk)
108///
109/// # Brief I/O
110///
111/// ```text
112/// VARIABLE I/O DESCRIPTION
113/// -------- --- --------------------------------------------------
114/// HANDLE I Handle assigned to the opened DSK file.
115/// CENTER I Central body ID code.
116/// SURFID I Surface ID code.
117/// DCLASS I Data class.
118/// FRAME I Reference frame.
119/// CORSYS I Coordinate system code.
120/// CORPAR I Coordinate system parameters.
121/// MNCOR1 I Minimum value of first coordinate.
122/// MXCOR1 I Maximum value of first coordinate.
123/// MNCOR2 I Minimum value of second coordinate.
124/// MXCOR2 I Maximum value of second coordinate.
125/// MNCOR3 I Minimum value of third coordinate.
126/// MXCOR3 I Maximum value of third coordinate.
127/// FIRST I Coverage start time.
128/// LAST I Coverage stop time.
129/// NV I Number of vertices.
130/// VRTCES I Vertices.
131/// NP I Number of plates.
132/// PLATES I Plates.
133/// SPAIXD I Double precision component of spatial index.
134/// SPAIXI I Integer component of spatial index.
135/// ANGMRG P Angular round-off margin.
136/// GENCLS P General surface DSK class.
137/// SVFCLS P Single-valued function DSK class.
138/// NSYPAR P Maximum number of coordinate system parameters in
139/// a DSK descriptor.
140/// MAXCGR P Maximum DSK type 2 coarse voxel count.
141/// MAXPLT P Maximum DSK type 2 plate count.
142/// MAXVOX P Maximum DSK type 2 voxel count.
143/// MAXVRT P Maximum DSK type 2 vertex count.
144/// ```
145///
146/// # Detailed Input
147///
148/// ```text
149/// HANDLE is the DAS file handle associated with a DSK file.
150/// The file must be open for write access.
151///
152/// CENTER is the ID code of the body whose surface is described
153/// by the input plate model. CENTER refers to an
154/// ephemeris object.
155///
156/// SURFID is the ID code of the surface described by the input
157/// plate model. Multiple surfaces (for example, surfaces
158/// having different resolutions) may be associated with
159/// a given body.
160///
161/// DCLASS is the data class of the input data set. See the
162/// INCLUDE file dskdsc.inc for values and meanings.
163///
164/// FRAME is the name of the reference frame with respect to
165/// which the input data are expressed.
166///
167/// CORSYS is the coordinate system in which the spatial coverage
168/// of the input data is expressed. CORSYS is an integer
169/// code. The supported values of CORPAR are
170///
171/// Parameter name Coordinate system
172/// -------------- -----------------
173/// LATSYS Planetocentric latitudinal
174/// RECSYS Rectangular (Cartesian)
175/// PDTSYS Planetodetic
176///
177/// See the INCLUDE file dskdsc.inc for parameter
178/// declarations.
179///
180///
181/// CORPAR is an array of parameters associated with the input
182/// coordinate system.
183///
184/// For latitudinal and rectangular coordinates, CORPAR
185/// is ignored.
186///
187/// For planetodetic coordinates, the contents of CORPAR
188/// are:
189///
190/// Element Contents
191/// --------- -----------------------------------
192/// CORPAR(1) Equatorial radius of reference
193/// spheroid.
194///
195/// CORPAR(2) Flattening coefficient. The polar
196/// radius of the reference spheroid
197/// is given by
198///
199/// CORPAR(1) * ( 1 - CORPAR(2) )
200///
201/// CORPAR(3)...
202/// CORPAR(NSYPAR) Unused.
203///
204///
205/// MNCOR1,
206/// MXCOR1,
207/// MNCOR2,
208/// MXCOR2,
209/// MNCOR3,
210/// MXCOR3 are, respectively, lower and upper bounds of
211/// each of the coordinates of the input data, where the
212/// coordinate system is defined by CORSYS and CORPAR.
213/// These bounds define the region for which the segment
214/// provides data.
215///
216/// Distance units are km. Angular units are radians.
217///
218/// The interpretation of these bounds depends on the data
219/// class; see DCLASS above.
220///
221/// Single-valued surfaces
222/// ----------------------
223///
224/// If the segment has data class SVFCLS (see
225/// dskdsc.inc), the segment defines a surface as a
226/// single-valued function of its domain coordinates:
227/// for example, it may define the radius of the
228/// surface as a function of planetocentric longitude
229/// and latitude, or Z as a function of X and Y.
230///
231/// In this case, the input data must cover a
232/// rectangle in dimensions 1 and 2 of the input
233/// coordinate system: the set of points
234///
235/// R = { (x,y): MNCOR1 <= x <= MXCOR1;
236/// MNCOR2 <= y <= MXCOR2 }
237///
238/// must be completely covered by the input data. In
239/// other words, for each point (x,y) of R, there must
240/// be some plate containing a point whose first two
241/// coordinates are (x,y).
242///
243/// The plate set may extend beyond the coordinate
244/// range defined by the bounds on the domain.
245///
246/// Normally for single-valued surfaces, MNCOR3 and
247/// MXCOR3 are the minimum and maximum values of the
248/// function attained on the domain.
249///
250///
251/// General surfaces
252/// ----------------
253///
254/// If the segment has data class GENCLS (see
255/// dskdsc.inc), the segment simply contains a
256/// collection of plates: no guarantees are made about
257/// the topology of the surface. The coordinate bounds
258/// simply indicate the spatial region for which the
259/// segment provides data.
260///
261/// Note that shapes of small bodies such as asteroids
262/// and comet nuclei may fall into the "general
263/// surface" category. Surface features such as cliffs,
264/// caves, and arches can prevent a surface from being
265/// represented as a single-valued function of latitude
266/// and longitude, for example.
267///
268///
269/// Longitude interpretation and restrictions
270/// -----------------------------------------
271///
272/// The following rules apply to longitudes provided in
273/// the arguments
274///
275/// MNCOR1
276/// MXCOR1
277///
278/// All angles have units of radians. The tolerance
279/// ANGMRG is used for the comparisons shown below.
280///
281/// 1) Longitudes must be in the range
282///
283/// -2*pi : 2*pi
284///
285/// Values that are out of range by no more than
286/// ANGMRG are bracketed to be in range.
287///
288///
289/// 2) It is acceptable for the longitude bounds to be
290/// out of order. If
291///
292/// MXCOR1 < MNCOR1
293///
294/// then either MXCOR1 is treated by the DSK
295/// subsystem as though it were MXCOR1 + 2*pi, or
296/// MNCOR1 is treated as MNCOR1 - 2*pi: whichever
297/// shift puts the bounds in the allowed range is
298/// made.
299///
300/// The input longitude bounds must not be equal.
301/// If the lower bound is greater than the upper
302/// bound, the difference between the bounds must
303/// not be an integer multiple of 2*pi.
304///
305/// 3) MXCOR1 must not exceed MNCOR1 by more than 2*pi.
306/// Values that are out of range by no more than
307/// ANGMRG are bracketed to be in range.
308///
309///
310/// FIRST,
311/// LAST are the endpoints of the time interval over which
312/// this data set is applicable. These times are
313/// expressed as seconds past J2000 TDB.
314///
315/// NV is the number of vertices belonging to the plate
316/// model.
317///
318/// VRTCES is an array of coordinates of the vertices.
319/// The Ith vertex occupies elements (1:3,I) of
320/// this array.
321///
322/// NP is the number of plates in the plate model.
323///
324/// PLATES is an array representing the plates of the model.
325/// The elements of PLATES are vertex indices. The vertex
326/// indices of the Ith plate occupy elements (1:3,I) of
327/// this array.
328///
329/// SPAIXD,
330/// SPAIXI are, respectively, the double precision and integer
331/// components of the spatial index of the segment.
332///
333/// It is strongly recommended that the helper routine
334/// DSKMI2 be used to create these arrays. See the
335/// $Examples section below.
336/// ```
337///
338/// # Detailed Output
339///
340/// ```text
341/// None. This routine operates by side effects.
342/// ```
343///
344/// # Parameters
345///
346/// ```text
347/// See the SPICELIB include files
348///
349/// dsk02.inc
350/// dskdsc.inc
351/// dsktol.inc
352///
353/// for declarations and detailed descriptions of the parameters
354/// referenced in this header.
355/// ```
356///
357/// # Exceptions
358///
359/// ```text
360/// 1) If the reference frame name FRAME could not be mapped to
361/// an ID code, the error SPICE(FRAMEIDNOTFOUND) is signaled.
362///
363/// 2) If the segment stop time precedes the start time, the
364/// error SPICE(TIMESOUTOFORDER) is signaled.
365///
366/// 3) If an input longitude value is outside the range
367///
368/// [ -2*pi - ANGMRG, 2*pi + ANGMRG ]
369///
370/// the error SPICE(VALUEOUTOFRANGE) is signaled. Longitudes
371/// outside of the range by a smaller amount than ANGMRG will be
372/// truncated to lie in the interval [-2*pi, 2*pi].
373///
374/// 4) If the absolute value of the difference between the input
375/// maximum longitude and the minimum longitude is more than 2*pi
376/// + ANGMRG, the error SPICE(INVALIDLONEXTENT) is signaled.
377/// If either longitude bound exceeds the other by an amount
378/// between 2*pi and 2*pi+ANGMRG, the larger value will be
379/// truncated to the smaller value plus 2*pi.
380///
381/// 5) If an input latitude value is outside the range
382///
383/// [ -pi/2 - ANGMRG, pi/2 + ANGMRG ]
384///
385/// the error SPICE(VALUEOUTOFRANGE) is signaled. Latitudes
386/// outside of the range by a smaller amount than ANGMRG will be
387/// truncated to lie in the interval [-pi/2, pi/2].
388///
389/// 6) If the coordinate system is latitudinal and the lower radius
390/// bound is negative, or if the upper radius bound is
391/// non-positive, the error SPICE(VALUEOUTOFRANGE) is signaled.
392///
393/// 7) If the coordinate system is latitudinal or planetodetic and
394/// the bounds of the latitude, radius or altitude coordinate are
395/// out of order, the error SPICE(BOUNDSOUTOFORDER) is signaled.
396///
397/// 8) If the coordinate system is latitudinal or planetodetic and
398/// the lower and upper bounds of the longitude, latitude, radius
399/// or altitude coordinate, respectively, are equal, the error
400/// SPICE(ZEROBOUNDSEXTENT) is signaled. If the lower
401/// longitude bound is greater than the upper bound, and if the
402/// difference between the bounds is an integer multiple of 2*pi,
403/// the same error is signaled.
404///
405/// 9) If the coordinate system is planetodetic and the input
406/// equatorial radius is non-positive, the error
407/// SPICE(VALUEOUTOFRANGE) is signaled.
408///
409/// 10) If the coordinate system is planetodetic and the input
410/// flattening coefficient is greater than or equal to 1, the
411/// error SPICE(VALUEOUTOFRANGE) is signaled.
412///
413/// 11) If the coordinate system is planetodetic, and if the minimum
414/// altitude is less than the maximum of
415///
416/// 2 2
417/// { -(B / A), -(A / B) }
418///
419/// where A and B are the semi-major and semi-minor axis lengths
420/// of the reference ellipsoid, the error
421/// SPICE(DEGENERATESURFACE) is signaled.
422///
423/// 12) If the coordinate system is rectangular and any coordinate
424/// lower bound is greater than or equal to the corresponding
425/// upper bound, the error SPICE(BOUNDSOUTOFORDER) is signaled.
426///
427/// 13) If the coordinate system code is not recognized, the error
428/// SPICE(NOTSUPPORTED) is signaled.
429///
430/// 14) If any vertex index belonging to an input plate is outside of
431/// the range 1:NV, the error SPICE(BADVERTEXINDEX) is signaled.
432///
433/// 15) If NV is less than 1 or greater than MAXVRT, the error
434/// SPICE(VALUEOUTOFRANGE) is signaled.
435///
436/// 16) If NP is less than 1 or greater than MAXPLT, the error
437/// SPICE(VALUEOUTOFRANGE) is signaled.
438///
439/// 17) If any voxel grid extent is less than 1 or greater than
440/// MAXVOX, the error SPICE(VALUEOUTOFRANGE) is signaled.
441///
442/// 18) If the voxel count is greater than MAXVOX, the error
443/// SPICE(VALUEOUTOFRANGE) is signaled.
444///
445/// 19) If the coarse voxel count is less than 1 or greater than
446/// MAXCGR, the error SPICE(VALUEOUTOFRANGE) is signaled.
447///
448/// 20) If the coarse voxel scale is less than 1 or more than
449/// the cube root of the fine voxel count, the error
450/// SPICE(VALUEOUTOFRANGE) is signaled.
451///
452/// 21) If the cube of the coarse voxel scale does not divide the
453/// fine voxel count evenly, the error SPICE(INCOMPATIBLESCALE)
454/// is signaled.
455///
456/// 22) If the input data class is not recognized, the error
457/// SPICE(NOTSUPPORTED) is signaled.
458///
459/// 23) If an error occurs while writing the segment to the output
460/// DSK file, the error is signaled by a routine in the call
461/// tree of this routine.
462/// ```
463///
464/// # Files
465///
466/// ```text
467/// See argument HANDLE.
468/// ```
469///
470/// # Particulars
471///
472/// ```text
473/// This routine writes a type 2 segment to a DSK file that
474/// has been opened for write access.
475///
476/// Users planning to create DSK files should consider whether the
477/// SPICE DSK creation utility MKDSK may be suitable for their needs.
478///
479/// This routine is supported by the routines DSKMI2 and DSKRB2
480/// DSKMI2 simplifies use of this routine by creating the "spatial
481/// index" arrays required as inputs by this routine. DSKRB2 computes
482/// bounds on the third coordinate of the input plate set.
483///
484/// Spatial Indexes
485/// ===============
486///
487/// A spatial index is a group of data structures that facilitates
488/// rapid high-level computations involving sets of plates. The data
489/// structures created by this routine are aggregated into arrays
490/// of type INTEGER and type DOUBLE PRECISION.
491///
492///
493/// Voxel grids
494/// ===========
495///
496/// A key geometric computation---probably the most important, as it
497/// serves as a foundation for other high-level computations---is
498/// finding the intersection of a ray with the plate set. DSK type 2
499/// segments use data structures called "voxel grids" as part of
500/// their indexing mechanism. There is a "coarse grid": a box that
501/// completely encloses a DSK type 2 segment's plate set, and which
502/// is composed of identically-sized cubes called "coarse voxels."
503/// Each coarse voxel in composed of smaller cubes called "fine
504/// voxels." When the term "voxel" is used without qualification, it
505/// refers to fine voxels.
506///
507/// Type 2 DSK segments contain data structures that associate plates
508/// with the fine voxels intersected by those plates. These
509/// structures enable the type 2 DSK software to rapidly find plates
510/// in a given region of space.
511///
512/// Voxel scales
513/// ============
514///
515/// There are two voxel scales:
516///
517/// - The coarse voxel scale is the integer ratio of the
518/// edge length of a coarse voxel to the edge length of
519/// a fine voxel
520///
521/// - The fine voxel scale is the double precision ratio
522/// of the edge length of a fine voxel to the average
523/// extent of the plates in the input plate set. "Extents"
524/// of a plate are the absolute values of the differences
525/// between the respective maximum and minimum X, Y, and Z
526/// coordinates of the plate's vertices.
527///
528/// Voxel scales determine the resolution of the voxel grid.
529/// Voxel scales must be chosen to satisfy size constraints and
530/// provide reasonable plate lookup performance.
531///
532/// The following considerations apply to spatial indexes of
533/// type 2 DSK segments:
534///
535/// 1) The maximum number of coarse voxels is fixed at MAXCGR
536/// (declared in dsk02.inc).
537///
538/// 2) If there are too few fine voxels, the average number of
539/// plates per fine voxel will be very large. This largely
540/// negates the performance improvement afforded by having an
541/// index. Also, the number of plates per voxel may exceed
542/// limits imposed by DSK subroutines that use static arrays.
543///
544/// 3) If there are too many fine voxels, the average number of
545/// voxels intersected by a given plate may be too large for
546/// all the plate-voxel associations to be stored. In
547/// addition, the time needed to examine the plate lists for
548/// each voxel (including the empty ones) may become quite
549/// large, again negating the value of the index.
550///
551/// In many cases, voxel scales yielding optimum performance must be
552/// determined by experiment. However, the following heuristics can
553/// provide reasonable starting values:
554///
555/// Let NP be the number of plates. Let FS be the fine voxel
556/// scale. Then a reasonable value of FS may be
557///
558/// (0.25D0)
559/// FS = NP / 8.D0
560///
561/// In general, FS should not smaller than 1.
562/// ```
563///
564/// # Examples
565///
566/// ```text
567/// The numerical results shown for this example may differ across
568/// platforms. The results depend on the SPICE kernels used as
569/// input, the compiler and supporting libraries, and the machine
570/// specific arithmetic implementation.
571///
572/// 1) Create a three-segment DSK file using plate model data for
573/// Phobos. Use latitudinal, rectangular, and planetodetic
574/// coordinates in the respective segments. This is not a
575/// realistic example, but it serves to demonstrate use of
576/// the supported coordinate systems.
577///
578/// Use the DSK kernel below to provide, for simplicity, the
579/// input plate and vertex data. The selected input file has one
580/// segment.
581///
582/// phobos_3_3.bds
583///
584///
585/// Example code begins here.
586///
587///
588/// C
589/// C Example program for DSKW02, DSKMI2, and DSKRB2
590/// C
591/// C Create a three-segment DSK file using plate model
592/// C data for Phobos. Use latitudinal, rectangular, and
593/// C planetodetic coordinates in the respective segments.
594/// C
595/// C For simplicity, use an existing DSK file to provide
596/// C the input plate and vertex data. The selected input
597/// C file has one segment.
598/// C
599/// C Version 1.0.0 22-JAN-2016 (NJB)
600/// C
601/// PROGRAM DSKW02_EX1
602/// IMPLICIT NONE
603///
604/// INCLUDE 'dla.inc'
605/// INCLUDE 'dskdsc.inc'
606/// INCLUDE 'dsk02.inc'
607///
608/// C
609/// C SPICELIB functions
610/// C
611/// DOUBLE PRECISION JYEAR
612/// DOUBLE PRECISION PI
613/// C
614/// C Local parameters
615/// C
616/// INTEGER FRNMLN
617/// PARAMETER ( FRNMLN = 32 )
618///
619/// INTEGER NSEG
620/// PARAMETER ( NSEG = 3 )
621///
622/// INTEGER NAMLEN
623/// PARAMETER ( NAMLEN = 20 )
624///
625/// INTEGER FILSIZ
626/// PARAMETER ( FILSIZ = 255 )
627///
628/// INTEGER LNSIZE
629/// PARAMETER ( LNSIZE = 80 )
630///
631/// INTEGER NCOR
632/// PARAMETER ( NCOR = 4 )
633///
634/// C
635/// C Local variables
636/// C
637/// CHARACTER*(NAMLEN) CORNAM ( NCOR )
638/// CHARACTER*(FILSIZ) DSK
639/// CHARACTER*(FRNMLN) FRAME
640/// CHARACTER*(FILSIZ) INDSK
641/// CHARACTER*(LNSIZE) LINE
642/// C
643/// C Note: the values of MAXVRT and MAXPLT declared
644/// C in dsk02.inc, and the integer spatial index
645/// C dimension SPAISZ are very large. Smaller buffers
646/// C can be used for most applications.
647/// C
648/// DOUBLE PRECISION CORPAR ( NSYPAR )
649/// DOUBLE PRECISION F
650/// DOUBLE PRECISION FINSCL
651/// DOUBLE PRECISION FIRST
652/// DOUBLE PRECISION LAST
653/// DOUBLE PRECISION MNCOR1
654/// DOUBLE PRECISION MNCOR2
655/// DOUBLE PRECISION MNCOR3
656/// DOUBLE PRECISION MXCOR1
657/// DOUBLE PRECISION MXCOR2
658/// DOUBLE PRECISION MXCOR3
659/// DOUBLE PRECISION RE
660/// DOUBLE PRECISION RP
661/// DOUBLE PRECISION SPAIXD ( IXDFIX )
662/// DOUBLE PRECISION VRTCES ( 3, MAXVRT )
663///
664/// INTEGER CENTER
665/// INTEGER CORSCL
666/// INTEGER CORSYS
667/// INTEGER DCLASS
668/// INTEGER DLADSC ( DLADSZ )
669/// INTEGER HANDLE
670/// INTEGER INHAN
671/// INTEGER NP
672/// INTEGER NV
673/// INTEGER PLATES ( 3, MAXPLT )
674/// INTEGER SEGNO
675/// INTEGER SPAIXI ( SPAISZ )
676/// INTEGER SURFID
677/// INTEGER VOXPSZ
678/// INTEGER VOXLSZ
679/// INTEGER WORK ( 2, MAXCEL )
680/// INTEGER WORKSZ
681///
682/// LOGICAL FOUND
683/// C
684/// C Saved variables
685/// C
686/// C Save all large arrays to avoid stack problems.
687/// C
688/// SAVE
689/// C
690/// C Initial values
691/// C
692/// DATA CORNAM / 'radius',
693/// . 'Z-coordinate',
694/// . 'Z-coordinate',
695/// . 'altitude' /
696///
697/// C
698/// C Assign names of input and output DSK files.
699/// C
700/// INDSK = 'phobos_3_3.bds'
701/// DSK = 'phobos_3_3_3seg.bds'
702/// C
703/// C Open input DSK for read access; find first segment.
704/// C
705/// CALL DASOPR ( INDSK, INHAN )
706/// CALL DLABFS ( INHAN, DLADSC, FOUND )
707/// C
708/// C Fetch vertices and plates from input DSK file.
709/// C
710/// WRITE (*,*) 'Reading input data...'
711///
712/// CALL DSKV02 ( INHAN, DLADSC, 1, MAXVRT, NV, VRTCES )
713/// CALL DSKP02 ( INHAN, DLADSC, 1, MAXPLT, NP, PLATES )
714///
715/// WRITE (*,*) 'Done.'
716/// C
717/// C Set input array sizes required by DSKMI2.
718/// C
719/// VOXPSZ = MAXVXP
720/// VOXLSZ = MXNVLS
721/// WORKSZ = MAXCEL
722/// C
723/// C Set fine and coarse voxel scales. (These usually
724/// C need to determined by experimentation.)
725/// C
726/// FINSCL = 5.D0
727/// CORSCL = 4
728/// C
729/// C Open a new DSK file.
730/// C
731/// CALL DSKOPN ( DSK, DSK, 0, HANDLE )
732/// C
733/// C Create three segments and add them to the file.
734/// C
735/// DO SEGNO = 1, NSEG
736/// C
737/// C Create spatial index.
738/// C
739/// WRITE (*,*) 'Creating segment ', SEGNO
740/// WRITE (*,*) 'Creating spatial index...'
741///
742/// CALL DSKMI2 ( NV, VRTCES, NP, PLATES, FINSCL,
743/// . CORSCL, WORKSZ, VOXPSZ, VOXLSZ, .TRUE.,
744/// . SPAISZ, WORK, SPAIXD, SPAIXI )
745///
746/// WRITE (*,*) 'Done.'
747/// C
748/// C Set up inputs describing segment attributes:
749/// C
750/// C - Central body: Phobos
751/// C - Surface ID code: user's choice.
752/// C We use the segment number here.
753/// C - Data class: general (arbitrary) shape
754/// C - Body-fixed reference frame
755/// C - Time coverage bounds (TBD)
756/// C
757/// CENTER = 401
758/// SURFID = SEGNO
759/// DCLASS = GENCLS
760/// FRAME = 'IAU_PHOBOS'
761///
762/// FIRST = -50 * JYEAR()
763/// LAST = 50 * JYEAR()
764/// C
765/// C Set the coordinate system and coordinate system
766/// C bounds based on the segment index.
767/// C
768/// C Zero out the coordinate parameters to start.
769/// C
770/// CALL CLEARD ( NSYPAR, CORPAR )
771///
772/// IF ( SEGNO .EQ. 1 ) THEN
773/// C
774/// C Use planetocentric latitudinal coordinates. Set
775/// C the longitude and latitude bounds.
776/// C
777/// CORSYS = LATSYS
778///
779/// MNCOR1 = -PI()
780/// MXCOR1 = PI()
781/// MNCOR2 = -PI()/2
782/// MXCOR2 = PI()/2
783///
784/// ELSE IF ( SEGNO .EQ. 2 ) THEN
785/// C
786/// C Use rectangular coordinates. Set the
787/// C X and Y bounds.
788/// C
789/// C The bounds shown here were derived from
790/// C the plate data. They lie slightly outside
791/// C of the range spanned by the plates.
792/// C
793/// CORSYS = RECSYS
794///
795/// MNCOR1 = -1.3D0
796/// MXCOR1 = 1.31D0
797/// MNCOR2 = -1.21D0
798/// MXCOR2 = 1.2D0
799///
800/// ELSE
801/// C
802/// C Set the coordinate system to planetodetic.
803/// C
804/// CORSYS = PDTSYS
805///
806/// MNCOR1 = -PI()
807/// MXCOR1 = PI()
808/// MNCOR2 = -PI()/2
809/// MXCOR2 = PI()/2
810/// C
811/// C We'll use equatorial and polar radii from
812/// C pck00010.tpc. These normally would be fetched
813/// C at run time, but for simplicity, we'll use
814/// C hard-coded values.
815///
816/// RE = 13.0D0
817/// RP = 9.1D0
818/// F = ( RE - RP ) / RE
819///
820/// CORPAR(1) = RE
821/// CORPAR(2) = F
822///
823/// END IF
824/// C
825/// C Compute plate model radius bounds.
826/// C
827/// LINE = 'Computing # bounds of plate set...'
828///
829/// CALL REPMC ( LINE, '#', CORNAM(CORSYS), LINE )
830/// WRITE (*,*) LINE
831///
832/// CALL DSKRB2 ( NV, VRTCES, NP, PLATES,
833/// . CORSYS, CORPAR, MNCOR3, MXCOR3 )
834///
835/// WRITE (*,*) 'Done.'
836/// C
837/// C Write the segment to the file.
838/// C
839/// WRITE (*,*) 'Writing segment...'
840///
841/// CALL DSKW02 ( HANDLE,
842/// . CENTER, SURFID, DCLASS, FRAME, CORSYS,
843/// . CORPAR, MNCOR1, MXCOR1, MNCOR2, MXCOR2,
844/// . MNCOR3, MXCOR3, FIRST, LAST, NV,
845/// . VRTCES, NP, PLATES, SPAIXD, SPAIXI )
846///
847/// WRITE (*,*) 'Done.'
848///
849/// END DO
850/// C
851/// C Segregate the data records in the DSK file and
852/// C close the file.
853/// C
854/// WRITE (*,*) 'Segregating and closing DSK file...'
855///
856/// CALL DSKCLS ( HANDLE, .TRUE. )
857///
858/// WRITE (*,*) 'Done.'
859/// END
860///
861///
862/// When this program was executed on a Mac/Intel/gfortran/64-bit
863/// platform, the output was:
864///
865///
866/// Reading input data...
867/// Done.
868/// Creating segment 1
869/// Creating spatial index...
870/// Done.
871/// Computing radius bounds of plate set...
872/// Done.
873/// Writing segment...
874/// Done.
875/// Creating segment 2
876/// Creating spatial index...
877/// Done.
878/// Computing Z-coordinate bounds of plate set...
879/// Done.
880/// Writing segment...
881/// Done.
882/// Creating segment 3
883/// Creating spatial index...
884/// Done.
885/// Computing altitude bounds of plate set...
886/// Done.
887/// Writing segment...
888/// Done.
889/// Segregating and closing DSK file...
890/// Done.
891///
892///
893/// Note that after run completion, a new DSK exists in the output
894/// directory.
895/// ```
896///
897/// # Author and Institution
898///
899/// ```text
900/// N.J. Bachman (JPL)
901/// J. Diaz del Rio (ODC Space)
902/// ```
903///
904/// # Version
905///
906/// ```text
907/// - SPICELIB Version 1.0.1, 28-AUG-2021 (JDR) (NJB)
908///
909/// Edited the header to comply with NAIF standard.
910/// Added solution to code example.
911///
912/// Corrected description of coverage requirements for class 1
913/// segments.
914///
915/// Deleted paragraph saying that, except for changes made to
916/// move longitude values into range, the values are stored in
917/// segment "as is."
918///
919/// - SPICELIB Version 1.0.0, 04-MAR-2017 (NJB)
920///
921/// Fixed some comment typos.
922///
923/// 10-OCT-2016 (NJB)
924///
925/// New error checks on inputs were added.
926///
927/// 07-MAR-2016 (NJB)
928///
929/// New error checks on inputs were added.
930///
931/// Argument list change: spatial index is now passed in
932/// as two arrays: SPAIXD and SPAIXI.
933///
934/// Argument CORPAR was added.
935///
936/// Double precision data are now written to the output
937/// segment before integer data.
938///
939/// 22-AUG-2012 (NJB)
940///
941/// Bug fix: corrected upper bound in test for
942/// vertex count.
943///
944/// 13-MAY-2010 (NJB)
945///
946/// Updated to reflect new type 2 segment design: normal
947/// vectors, plate centers, and lengths of longest plate sides
948/// are no longer stored in these segments.
949///
950/// 03-APR-2010 (NJB)
951///
952/// New interface; general coordinates are supported. Time
953/// bounds, surface ID, data class, and bounds of third
954/// coordinate have been added. Albedo inputs have been
955/// deleted.
956///
957/// 09-OCT-2009 (NJB)
958///
959/// Header was added.
960///
961/// 31-OCT-2006 (NJB)
962///
963/// Input arguments CGSCAL and VTXBDS were removed.
964///
965/// 27-JUN-2006 (NJB)
966///
967/// Initial version.
968/// ```
969pub fn dskw02(
970 ctx: &mut SpiceContext,
971 handle: i32,
972 center: i32,
973 surfid: i32,
974 dclass: i32,
975 frame: &str,
976 corsys: i32,
977 corpar: &[f64],
978 mncor1: f64,
979 mxcor1: f64,
980 mncor2: f64,
981 mxcor2: f64,
982 mncor3: f64,
983 mxcor3: f64,
984 first: f64,
985 last: f64,
986 nv: i32,
987 vrtces: &[[f64; 3]],
988 np: i32,
989 plates: &[[i32; 3]],
990 spaixd: &[f64],
991 spaixi: &[i32],
992) -> crate::Result<()> {
993 DSKW02(
994 handle,
995 center,
996 surfid,
997 dclass,
998 frame.as_bytes(),
999 corsys,
1000 corpar,
1001 mncor1,
1002 mxcor1,
1003 mncor2,
1004 mxcor2,
1005 mncor3,
1006 mxcor3,
1007 first,
1008 last,
1009 nv,
1010 vrtces.as_flattened(),
1011 np,
1012 plates.as_flattened(),
1013 spaixd,
1014 spaixi,
1015 ctx.raw_context(),
1016 )?;
1017 ctx.handle_errors()?;
1018 Ok(())
1019}
1020
1021//$Procedure DSKW02 ( DSK, write type 2 segment )
1022pub fn DSKW02(
1023 HANDLE: i32,
1024 CENTER: i32,
1025 SURFID: i32,
1026 DCLASS: i32,
1027 FRAME: &[u8],
1028 CORSYS: i32,
1029 CORPAR: &[f64],
1030 MNCOR1: f64,
1031 MXCOR1: f64,
1032 MNCOR2: f64,
1033 MXCOR2: f64,
1034 MNCOR3: f64,
1035 MXCOR3: f64,
1036 FIRST: f64,
1037 LAST: f64,
1038 NV: i32,
1039 VRTCES: &[f64],
1040 NP: i32,
1041 PLATES: &[i32],
1042 SPAIXD: &[f64],
1043 SPAIXI: &[i32],
1044 ctx: &mut Context,
1045) -> f2rust_std::Result<()> {
1046 let CORPAR = DummyArray::new(CORPAR, 1..);
1047 let VRTCES = DummyArray2D::new(VRTCES, 1..=3, 1..=NV);
1048 let PLATES = DummyArray2D::new(PLATES, 1..=3, 1..=NP);
1049 let SPAIXD = DummyArray::new(SPAIXD, 1..);
1050 let SPAIXI = DummyArray::new(SPAIXI, 1..);
1051 let mut A: f64 = 0.0;
1052 let mut ALTLIM: f64 = 0.0;
1053 let mut B: f64 = 0.0;
1054 let mut DESCR = StackArray::<f64, 24>::new(1..=DSKDSZ);
1055 let mut R: f64 = 0.0;
1056 let mut SEGBDS = StackArray2D::<f64, 4>::new(1..=2, 1..=2);
1057 let mut VOXORI = StackArray::<f64, 3>::new(1..=3);
1058 let mut VOXSIZ: f64 = 0.0;
1059 let mut VTXBDS = StackArray2D::<f64, 6>::new(1..=2, 1..=3);
1060 let mut CGRSCL: i32 = 0;
1061 let mut FRMCDE: i32 = 0;
1062 let mut K: i32 = 0;
1063 let mut NCGR: i32 = 0;
1064 let mut NVXTOT: i32 = 0;
1065 let mut PVOXPL: i32 = 0;
1066 let mut PVOXPT: i32 = 0;
1067 let mut PVTXPL: i32 = 0;
1068 let mut PVTXPT: i32 = 0;
1069 let mut Q: i32 = 0;
1070 let mut VGREXT = StackArray::<i32, 3>::new(1..=3);
1071 let mut VOXNPT: i32 = 0;
1072 let mut VOXNPL: i32 = 0;
1073 let mut VTXNPL: i32 = 0;
1074
1075 //
1076 // SPICELIB functions
1077 //
1078
1079 //
1080 // Local parameters
1081 //
1082
1083 //
1084 // Local variables
1085 //
1086
1087 if RETURN(ctx) {
1088 return Ok(());
1089 }
1090
1091 CHKIN(b"DSKW02", ctx)?;
1092
1093 //
1094 // Map the input reference frame name to an ID code.
1095 //
1096 NAMFRM(FRAME, &mut FRMCDE, ctx)?;
1097
1098 if (FRMCDE == 0) {
1099 SETMSG(b"Input reference frame # could not be mapped to an ID code. The frame name might be misspelled, or possibly a required frame kernel was not loaded. ", ctx);
1100 ERRCH(b"#", FRAME, ctx);
1101 SIGERR(b"SPICE(FRAMEIDNOTFOUND)", ctx)?;
1102 CHKOUT(b"DSKW02", ctx)?;
1103 return Ok(());
1104 }
1105
1106 //
1107 // Make sure the time bounds are in order.
1108 //
1109 if (LAST <= FIRST) {
1110 SETMSG(
1111 b"Segment time bounds must be increasing; bounds were #:#.",
1112 ctx,
1113 );
1114 ERRDP(b"#", FIRST, ctx);
1115 ERRDP(b"#", LAST, ctx);
1116 SIGERR(b"SPICE(TIMESOUTOFORDER)", ctx)?;
1117 CHKOUT(b"DSKW02", ctx)?;
1118 return Ok(());
1119 }
1120
1121 //
1122 // If applicable, check segment boundaries. Check the
1123 // coordinate system as well.
1124 //
1125 if ((CORSYS == LATSYS) || (CORSYS == PDTSYS)) {
1126 //
1127 // Reject invalid latitudes and longitudes. Move
1128 // values that are slightly out of range into range.
1129 //
1130 // Longitude bounds must be distinct.
1131 //
1132 if (MNCOR1 == MXCOR1) {
1133 SETMSG(b"Minimum longitude # radians (# degrees) was equal to maximum longitude. Longitude bounds must be distinct.", ctx);
1134 ERRDP(b"#", MNCOR1, ctx);
1135 ERRDP(b"#", (MNCOR1 * DPR(ctx)), ctx);
1136 SIGERR(b"SPICE(ZEROBOUNDSEXTENT)", ctx)?;
1137 CHKOUT(b"DSKW02", ctx)?;
1138 return Ok(());
1139 }
1140
1141 //
1142 // Check minimum longitude.
1143 //
1144 if ((MNCOR1 < (-TWOPI(ctx) - ANGMRG)) || (MNCOR1 > (TWOPI(ctx) - ANGMRG))) {
1145 SETMSG(b"Minimum longitude # radians (# degrees) was outside of valid range [-2*pi, 2*pi - ANGMRG]", ctx);
1146 ERRDP(b"#", MNCOR1, ctx);
1147 ERRDP(b"#", (MNCOR1 * DPR(ctx)), ctx);
1148 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1149 CHKOUT(b"DSKW02", ctx)?;
1150 return Ok(());
1151 }
1152
1153 //
1154 // The minimum longitude is too small by ANGMRG, at worst.
1155 // Make it greater than or equal to -2*pi.
1156 //
1157 SEGBDS[[1, 1]] = intrinsics::DMAX1(&[-TWOPI(ctx), MNCOR1]);
1158
1159 //
1160 // Check maximum longitude.
1161 //
1162 if ((MXCOR1 < (-TWOPI(ctx) + ANGMRG)) || (MXCOR1 > (TWOPI(ctx) + ANGMRG))) {
1163 SETMSG(b"Maximum longitude # radians (# degrees) was outside of valid range [-2*pi+ANGMRG, 2*pi]", ctx);
1164 ERRDP(b"#", MXCOR1, ctx);
1165 ERRDP(b"#", (MXCOR1 * DPR(ctx)), ctx);
1166 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1167 CHKOUT(b"DSKW02", ctx)?;
1168 return Ok(());
1169 }
1170
1171 //
1172 // The maximum longitude is too large by ANGMRG, at worst.
1173 // Make it less than or equal to 2*pi.
1174 //
1175 SEGBDS[[2, 1]] = intrinsics::DMIN1(&[TWOPI(ctx), MXCOR1]);
1176
1177 //
1178 // The longitude extent cannot exceed 2*pi.
1179 //
1180 if (MXCOR1 > ((MNCOR1 + TWOPI(ctx)) + ANGMRG)) {
1181 SETMSG(
1182 b"Longitude bounds #:# radians (#:# degrees) are too far apart.",
1183 ctx,
1184 );
1185 ERRDP(b"#", MXCOR1, ctx);
1186 ERRDP(b"#", MXCOR2, ctx);
1187 ERRDP(b"#", (MXCOR1 * DPR(ctx)), ctx);
1188 ERRDP(b"#", (MXCOR2 * DPR(ctx)), ctx);
1189 SIGERR(b"SPICE(INVALIDLONEXTENT)", ctx)?;
1190 CHKOUT(b"DSKW02", ctx)?;
1191 return Ok(());
1192 }
1193
1194 if (MXCOR1 < ((MNCOR1 - TWOPI(ctx)) - ANGMRG)) {
1195 SETMSG(
1196 b"Longitude bounds #:# radians (#:# degrees) are too far apart.",
1197 ctx,
1198 );
1199 ERRDP(b"#", MXCOR1, ctx);
1200 ERRDP(b"#", MXCOR2, ctx);
1201 ERRDP(b"#", (MXCOR1 * DPR(ctx)), ctx);
1202 ERRDP(b"#", (MXCOR2 * DPR(ctx)), ctx);
1203 SIGERR(b"SPICE(INVALIDLONEXTENT)", ctx)?;
1204 CHKOUT(b"DSKW02", ctx)?;
1205 return Ok(());
1206 }
1207
1208 if (SEGBDS[[2, 1]] > SEGBDS[[1, 1]]) {
1209 //
1210 // The upper bound exceeds the lower by at most 2*pi + ANGMRG.
1211 // Trim the upper bound to make the difference no more than
1212 // 2*pi.
1213 //
1214 SEGBDS[[2, 1]] = intrinsics::DMIN1(&[SEGBDS[[2, 1]], (SEGBDS[[1, 1]] + TWOPI(ctx))]);
1215 } else if (SEGBDS[[2, 1]] < SEGBDS[[1, 1]]) {
1216 //
1217 // The lower bound exceeds the upper by at most 2*pi + ANGMRG.
1218 // Advance the upper bound, if necessary, to make the
1219 // difference no more than 2*pi.
1220 //
1221 SEGBDS[[2, 1]] = intrinsics::DMAX1(&[SEGBDS[[2, 1]], (SEGBDS[[1, 1]] - TWOPI(ctx))]);
1222 }
1223
1224 //
1225 // Make sure the adjusted longitude bounds don't describe an
1226 // interval that could be interpreted as having length zero,
1227 // if the bounds were placed in order. If the lower bound is
1228 // greater than the upper bound, then the difference between
1229 // the bounds must not be an integer multiple of 2*pi.
1230 //
1231 if ((SEGBDS[[2, 1]] == SEGBDS[[1, 1]]) || (SEGBDS[[2, 1]] == (SEGBDS[[1, 1]] - TWOPI(ctx))))
1232 {
1233 SETMSG(b"After adjustment, minimum longitude # radians (# degrees) was equal to maximum longitude. Longitude bounds must be distinct.", ctx);
1234 ERRDP(b"#", SEGBDS[[1, 1]], ctx);
1235 ERRDP(b"#", (MNCOR1 * DPR(ctx)), ctx);
1236 SIGERR(b"SPICE(ZEROBOUNDSEXTENT)", ctx)?;
1237 CHKOUT(b"DSKW02", ctx)?;
1238 return Ok(());
1239 }
1240
1241 //
1242 // Check minimum latitude.
1243 //
1244 if ((MNCOR2 < (-HALFPI(ctx) - ANGMRG)) || (MNCOR2 > (HALFPI(ctx) - ANGMRG))) {
1245 SETMSG(b"Minimum latitude # radians (# degrees) was outside of valid range [-pi/2, pi/2 - ANGMRG]", ctx);
1246 ERRDP(b"#", MNCOR2, ctx);
1247 ERRDP(b"#", (MNCOR2 * DPR(ctx)), ctx);
1248 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1249 CHKOUT(b"DSKW02", ctx)?;
1250 return Ok(());
1251 }
1252
1253 //
1254 // Trim the lower latitude bound to make it at least -pi/2.
1255 //
1256 SEGBDS[[1, 2]] = intrinsics::DMAX1(&[-HALFPI(ctx), MNCOR2]);
1257
1258 //
1259 // Check maximum latitude.
1260 //
1261 if ((MXCOR2 < (-HALFPI(ctx) + ANGMRG)) || (MXCOR2 > (HALFPI(ctx) + ANGMRG))) {
1262 SETMSG(b"Maximum latitude # radians (# degrees) was outside of valid range [-pi/2+ANGMRG, pi/2]", ctx);
1263 ERRDP(b"#", MXCOR2, ctx);
1264 ERRDP(b"#", (MXCOR2 * DPR(ctx)), ctx);
1265 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1266 CHKOUT(b"DSKW02", ctx)?;
1267 return Ok(());
1268 }
1269
1270 //
1271 // Trim the upper latitude bound to make it no more than -pi/2.
1272 //
1273 SEGBDS[[1, 2]] = intrinsics::DMAX1(&[-HALFPI(ctx), MNCOR2]);
1274 SEGBDS[[2, 2]] = intrinsics::DMIN1(&[HALFPI(ctx), MXCOR2]);
1275
1276 //
1277 // The latitude bounds must be in order.
1278 //
1279 if (MXCOR2 < MNCOR2) {
1280 SETMSG(b"Latitude bounds # and # are out of order.", ctx);
1281 ERRDP(b"#", MNCOR2, ctx);
1282 ERRDP(b"#", MXCOR2, ctx);
1283 SIGERR(b"SPICE(BOUNDSOUTOFORDER)", ctx)?;
1284 CHKOUT(b"DSKW02", ctx)?;
1285 return Ok(());
1286 }
1287
1288 if (CORSYS == LATSYS) {
1289 //
1290 // The coordinate system is latitudinal. Check radius
1291 // bounds.
1292 //
1293 if (MNCOR3 < 0.0) {
1294 SETMSG(b"Radius lower bound must be non-negative but was #.", ctx);
1295 ERRDP(b"#", MNCOR3, ctx);
1296 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1297 CHKOUT(b"DSKW02", ctx)?;
1298 return Ok(());
1299 }
1300
1301 if (MXCOR3 <= 0.0) {
1302 SETMSG(
1303 b"Radius upper bound must be strictly positive but was #.",
1304 ctx,
1305 );
1306 ERRDP(b"#", MXCOR3, ctx);
1307 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1308 CHKOUT(b"DSKW02", ctx)?;
1309 return Ok(());
1310 }
1311 }
1312
1313 if (CORSYS == PDTSYS) {
1314 //
1315 // The coordinate system is planetodetic. Check the coordinate
1316 // parameters as well.
1317 //
1318 if (CORPAR[1] <= 0.0) {
1319 SETMSG(
1320 b"Equatorial radius was #; this radius must be strictly positive.",
1321 ctx,
1322 );
1323 ERRDP(b"#", CORPAR[1], ctx);
1324 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1325 CHKOUT(b"DSKW02", ctx)?;
1326 return Ok(());
1327 }
1328
1329 if (CORPAR[2] >= 1.0) {
1330 SETMSG(
1331 b"Flattening coefficient was #; this value must be strictly less than 1.",
1332 ctx,
1333 );
1334 ERRDP(b"#", CORPAR[2], ctx);
1335 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1336 CHKOUT(b"DSKW02", ctx)?;
1337 return Ok(());
1338 }
1339
1340 //
1341 // Make sure the surface of minimum altitude is smooth and
1342 // non-self-intersecting.
1343 //
1344 A = CORPAR[1];
1345 B = (A * ((1 as f64) - CORPAR[2]));
1346
1347 ALTLIM = intrinsics::DMAX1(&[-(f64::powi(A, 2) / B), -(f64::powi(B, 2) / A)]);
1348
1349 if (MNCOR3 <= ALTLIM) {
1350 SETMSG(b"Reference ellipsoid has semi-axis lengths # and #. The minimum altitude was #. The minimum altitude is required to be greater than the maximum of {-(A**2)/B, -(B**2)/A}, which is #.", ctx);
1351 ERRDP(b"#", A, ctx);
1352 ERRDP(b"#", B, ctx);
1353 ERRDP(b"#", MNCOR3, ctx);
1354 ERRDP(b"#", ALTLIM, ctx);
1355 SIGERR(b"SPICE(DEGENERATESURFACE)", ctx)?;
1356 CHKOUT(b"DSKW02", ctx)?;
1357 return Ok(());
1358 }
1359 }
1360
1361 //
1362 // The bounds of the third coordinate, whether radius or altitude,
1363 // must be in order and must have positive extent.
1364 //
1365 if (MXCOR3 < MNCOR3) {
1366 if (CORSYS == LATSYS) {
1367 SETMSG(b"Radius bounds # and # are out of order", ctx);
1368 } else {
1369 SETMSG(b"Altitude bounds # and # are out of order.", ctx);
1370 }
1371
1372 ERRDP(b"#", MNCOR3, ctx);
1373 ERRDP(b"#", MXCOR3, ctx);
1374 SIGERR(b"SPICE(BOUNDSOUTOFORDER)", ctx)?;
1375 CHKOUT(b"DSKW02", ctx)?;
1376 return Ok(());
1377 }
1378
1379 if (MXCOR3 == MNCOR3) {
1380 if (CORSYS == LATSYS) {
1381 SETMSG(
1382 b"Radius bounds # and # must have positive extent but are equal.",
1383 ctx,
1384 );
1385 } else {
1386 SETMSG(
1387 b"Altitude bounds # and # must have positive extent but are equal.",
1388 ctx,
1389 );
1390 }
1391
1392 ERRDP(b"#", MNCOR3, ctx);
1393 ERRDP(b"#", MXCOR3, ctx);
1394 SIGERR(b"SPICE(ZEROBOUNDSEXTENT)", ctx)?;
1395 CHKOUT(b"DSKW02", ctx)?;
1396 return Ok(());
1397 }
1398 } else if (CORSYS == RECSYS) {
1399 //
1400 // All coordinate bounds must be in strictly increasing order.
1401 //
1402 if (((MXCOR1 <= MNCOR1) || (MXCOR2 <= MNCOR2)) || (MXCOR3 <= MNCOR3)) {
1403 SETMSG(b"Rectangular coordinate bounds must be strictly increasing in each dimension. The bounds were: X = #:#; Y = #:#; Z = #:#.", ctx);
1404 ERRDP(b"#", MNCOR1, ctx);
1405 ERRDP(b"#", MXCOR1, ctx);
1406 ERRDP(b"#", MNCOR2, ctx);
1407 ERRDP(b"#", MXCOR2, ctx);
1408 ERRDP(b"#", MNCOR3, ctx);
1409 ERRDP(b"#", MXCOR3, ctx);
1410 SIGERR(b"SPICE(BOUNDSOUTOFORDER)", ctx)?;
1411 CHKOUT(b"DSKW02", ctx)?;
1412 return Ok(());
1413 }
1414
1415 SEGBDS[[1, 1]] = MNCOR1;
1416 SEGBDS[[2, 1]] = MXCOR1;
1417 SEGBDS[[1, 2]] = MNCOR2;
1418 SEGBDS[[2, 2]] = MXCOR2;
1419 } else {
1420 SETMSG(b"Coordinate system code # is not recognized.", ctx);
1421 ERRINT(b"#", CORSYS, ctx);
1422 SIGERR(b"SPICE(NOTSUPPORTED)", ctx)?;
1423 CHKOUT(b"DSKW02", ctx)?;
1424 return Ok(());
1425 }
1426
1427 //
1428 // Check the data class.
1429 //
1430 if ((DCLASS < 1) || (DCLASS > 2)) {
1431 SETMSG(b"Data class # is not recognized.", ctx);
1432 ERRINT(b"#", DCLASS, ctx);
1433 SIGERR(b"SPICE(NOTSUPPORTED)", ctx)?;
1434 CHKOUT(b"DSKW02", ctx)?;
1435 return Ok(());
1436 }
1437
1438 //
1439 // Check NV and NP.
1440 //
1441 // Note that we don't apply Euler's law, since the data
1442 // set need not represent a complete surface.
1443 //
1444 if ((NV < 1) || (NV > MAXVRT)) {
1445 SETMSG(b"Vertex count NV = #; count must be in the range 1:#.", ctx);
1446 ERRINT(b"#", NV, ctx);
1447 ERRINT(b"#", MAXVRT, ctx);
1448 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1449 CHKOUT(b"DSKW02", ctx)?;
1450 return Ok(());
1451 }
1452
1453 if ((NP < 1) || (NP > MAXPLT)) {
1454 SETMSG(b"Plate count NP = #; count must be in the range 1:#.", ctx);
1455 ERRINT(b"#", NP, ctx);
1456 ERRINT(b"#", MAXPLT, ctx);
1457 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1458 CHKOUT(b"DSKW02", ctx)?;
1459 return Ok(());
1460 }
1461
1462 //
1463 // Check the vertex indices in the plates.
1464 //
1465 for I in 1..=NP {
1466 for J in 1..=3 {
1467 K = PLATES[[J, I]];
1468
1469 if ((K < 1) || (K > NV)) {
1470 SETMSG(b"Vertex index # of plate # was #; vertex indices must be in the range 1:NV. The input NV = #.", ctx);
1471 ERRINT(b"#", J, ctx);
1472 ERRINT(b"#", I, ctx);
1473 ERRINT(b"#", K, ctx);
1474 ERRINT(b"#", NV, ctx);
1475 SIGERR(b"SPICE(BADVERTEXINDEX)", ctx)?;
1476 CHKOUT(b"DSKW02", ctx)?;
1477 return Ok(());
1478 }
1479 }
1480 }
1481
1482 //
1483 // Locate the spatial index elements. Some of the elements are at
1484 // fixed addresses; for others the addresses must be calculated.
1485 //
1486 // The two components of the spatial index together contain the
1487 // following items:
1488 //
1489 // VGREXT is an array containing the extents of the voxel
1490 // grid in the X, Y, and Z directions of the
1491 // body-fixed frame. The extents are measured as
1492 // voxel counts.
1493 //
1494 // ORIGIN is the position, in the body-fixed, body-centered
1495 // reference frame associated with BODY, of the
1496 // origin of the both the fine and coarse voxel grids
1497 // associated with this model.
1498 //
1499 // VOXSIZ is the voxel edge length in km.
1500 //
1501 // CGRSCL is the coarse voxel grid scale factor: the edge
1502 // length of each coarse voxel is scaled up from the
1503 // length of a fine voxel edge by this factor.
1504 //
1505 // CGRPTR is an array of pointers associated with this
1506 // model's coarse voxel grid; these pointers map
1507 // one-dimensional coarse voxel indices to start
1508 // indices in the fine voxel pointer list.
1509 //
1510 // VOXNPT is the cardinality of the fine voxel pointer list.
1511 //
1512 // VOXPTR is the fine voxel pointer list. For each fine
1513 // voxel belonging to a non-empty coarse voxel, there
1514 // is a pointer in this list that identifies the
1515 // start index in VOXPLT of the list of plate indices
1516 // associated with this fine voxel.
1517 //
1518 // The start index in VOXPTR of the set of pointers
1519 // associated with a coarse voxel is given by the
1520 // element of CGRPTR associated with that coarse
1521 // voxel.
1522 //
1523 // Within a given coarse voxel, each fine voxel has
1524 // an associated one-dimensional offset from the
1525 // corner of the coarse voxel closest to the origin
1526 // of the voxel grids. This offset gives the location
1527 // in VOXPTR of the plate list pointer for the fine
1528 // voxel.
1529 //
1530 // VOXNPL is the cardinality of the plate list of the fine
1531 // voxel-plate mapping.
1532 //
1533 // VOXPLT is the plate list of the fine voxel-plate mapping.
1534 //
1535 // VTXPTR is the vertex pointer list.
1536 //
1537 // VTXNPL is the cardinality of the plate list of the
1538 // vertex-plate mapping.
1539 //
1540 //
1541 //
1542 // Extract double precision elements of the spatial index.
1543 //
1544 MOVED(SPAIXD.subarray(SIVTBD), 6, VTXBDS.as_slice_mut());
1545 VEQU(SPAIXD.subarray(SIVXOR), VOXORI.as_slice_mut());
1546
1547 VOXSIZ = SPAIXD[SIVXSZ];
1548
1549 //
1550 // Extract scalars and small fixed-size arrays from the integer
1551 // component of the spatial index.
1552 //
1553 // Fetch grid extents (in units of whole voxels):
1554 //
1555 MOVEI(SPAIXI.subarray(SIVGRX), 3, VGREXT.as_slice_mut());
1556
1557 //
1558 // Fetch coarse grid scale, voxel pointer count, and voxel-plate
1559 // list count.
1560 //
1561 CGRSCL = SPAIXI[SICGSC];
1562 VOXNPT = SPAIXI[SIVXNP];
1563 VOXNPL = SPAIXI[SIVXNL];
1564
1565 //
1566 // Create a pointer to the voxel-plate pointer array.
1567 //
1568 PVOXPT = (SICGRD + MAXCGR);
1569
1570 //
1571 // Create a pointer to the voxel-plate list array.
1572 //
1573 PVOXPL = (PVOXPT + VOXNPT);
1574
1575 //
1576 // Create a pointer to the vertex pointer array.
1577 //
1578 PVTXPT = (PVOXPL + VOXNPL);
1579
1580 //
1581 // Create a pointer to the vertex-plate list array.
1582 //
1583 PVTXPL = (PVTXPT + NV);
1584
1585 //
1586 // Fetch vertex-plate list size.
1587 //
1588 VTXNPL = SPAIXI[SIVTNL];
1589
1590 //
1591 // Check the input parameters.
1592 //
1593 //
1594 //
1595 // Make sure the voxel grid extents are within range.
1596 //
1597 for I in 1..=3 {
1598 if ((VGREXT[I] < 1) || (VGREXT[I] > MAXVOX)) {
1599 SETMSG(
1600 b"Voxel grid extents are = (#, #, #); all be in the range 1:#.",
1601 ctx,
1602 );
1603 ERRINT(b"#", VGREXT[1], ctx);
1604 ERRINT(b"#", VGREXT[2], ctx);
1605 ERRINT(b"#", VGREXT[3], ctx);
1606 ERRINT(b"#", MAXVOX, ctx);
1607 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1608 CHKOUT(b"DSKW02", ctx)?;
1609 return Ok(());
1610 }
1611 }
1612
1613 //
1614 // Make sure the number of voxels NVXTOT is within range.
1615 //
1616 NVXTOT = ((VGREXT[1] * VGREXT[2]) * VGREXT[3]);
1617
1618 if (NVXTOT > MAXVOX) {
1619 SETMSG(
1620 b"Fine voxel count NVXTOT = #; count must be in the range 1:#.",
1621 ctx,
1622 );
1623 ERRINT(b"#", NVXTOT, ctx);
1624 ERRINT(b"#", MAXVOX, ctx);
1625 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1626 CHKOUT(b"DSKW02", ctx)?;
1627 return Ok(());
1628 }
1629
1630 //
1631 // Check the coarse voxel scale. It must be at least 1, and its
1632 // cube must not exceed the fine voxel count.
1633 //
1634 if ((CGRSCL < 1) || ((CGRSCL as f64) > f64::powf(NVXTOT as f64, (1.0 / 3 as f64)))) {
1635 SETMSG(b"Coarse voxel scale = #; scale must be in the range 1:NVXTOT**3, where NVXTOT is the total fine voxel count. In this case, NVXTOT = #.", ctx);
1636 ERRINT(b"#", CGRSCL, ctx);
1637 ERRINT(b"#", NVXTOT, ctx);
1638 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1639 CHKOUT(b"DSKW02", ctx)?;
1640 return Ok(());
1641 }
1642
1643 //
1644 // The cube of the coarse scale must divide the total voxel count
1645 // evenly.
1646 //
1647 Q = (NVXTOT / intrinsics::pow(CGRSCL, 3));
1648 R = (NVXTOT - (Q * intrinsics::pow(CGRSCL, 3))) as f64;
1649
1650 if (R != 0 as f64) {
1651 SETMSG(b"Coarse voxel scale = #; the cube of the scale must divide NVXTOT evenly, where NVXTOT is the total fine voxel count. In this case, NVXTOT = #.", ctx);
1652 ERRINT(b"#", CGRSCL, ctx);
1653 ERRINT(b"#", NVXTOT, ctx);
1654 SIGERR(b"SPICE(INCOMPATIBLESCALE)", ctx)?;
1655 CHKOUT(b"DSKW02", ctx)?;
1656 return Ok(());
1657 }
1658
1659 //
1660 // NCGR is the number of voxels in the coarse voxel grid
1661 // associated with this model. Since each coarse voxel
1662 // is a cube containing an integer number of fine
1663 // voxels, this number is determined by NVXTOT and
1664 // CGRSCL.
1665 //
1666 NCGR = (NVXTOT / intrinsics::pow(CGRSCL, 3));
1667
1668 //
1669 // Make sure NCGR is within range.
1670 //
1671 if ((NCGR < 1) || (NCGR > MAXCGR)) {
1672 SETMSG(
1673 b"Coarse voxel count = #; count must be in the range 1:#.",
1674 ctx,
1675 );
1676 ERRINT(b"#", NCGR, ctx);
1677 ERRINT(b"#", MAXCGR, ctx);
1678 SIGERR(b"SPICE(VALUEOUTOFRANGE)", ctx)?;
1679 CHKOUT(b"DSKW02", ctx)?;
1680 return Ok(());
1681 }
1682
1683 //
1684 // Start a new DLA segment.
1685 //
1686 DLABNS(HANDLE, ctx)?;
1687
1688 if FAILED(ctx) {
1689 CHKOUT(b"DSKW02", ctx)?;
1690 return Ok(());
1691 }
1692
1693 //
1694 // Write the d.p. data to the segment first. In the segregated
1695 // segment, d.p. data will precede integer data, so less
1696 // rearrangement will occur this way.
1697 //
1698 //
1699 // First, fill in the DSK segment descriptor.
1700 //
1701 CLEARD(DSKDSZ, DESCR.as_slice_mut());
1702
1703 DESCR[SRFIDX] = SURFID as f64;
1704 DESCR[CTRIDX] = CENTER as f64;
1705 DESCR[CLSIDX] = DCLASS as f64;
1706 DESCR[TYPIDX] = 2 as f64;
1707 DESCR[FRMIDX] = FRMCDE as f64;
1708 DESCR[SYSIDX] = CORSYS as f64;
1709
1710 MOVED(CORPAR.as_slice(), NSYPAR, DESCR.subarray_mut(PARIDX));
1711
1712 DESCR[MN1IDX] = SEGBDS[[1, 1]];
1713 DESCR[MX1IDX] = SEGBDS[[2, 1]];
1714 DESCR[MN2IDX] = SEGBDS[[1, 2]];
1715 DESCR[MX2IDX] = SEGBDS[[2, 2]];
1716 DESCR[MN3IDX] = MNCOR3;
1717 DESCR[MX3IDX] = MXCOR3;
1718 DESCR[BTMIDX] = FIRST;
1719 DESCR[ETMIDX] = LAST;
1720
1721 //
1722 // Now write the descriptor into the segment.
1723 //
1724 DASADD(HANDLE, DSKDSZ, DESCR.as_slice(), ctx)?;
1725
1726 //
1727 // Add the voxel grid origin and voxel size.
1728 // Finish with the vertex data.
1729 //
1730 DASADD(HANDLE, 6, VTXBDS.as_slice(), ctx)?;
1731 DASADD(HANDLE, 3, VOXORI.as_slice(), ctx)?;
1732 DASADD(HANDLE, 1, &[VOXSIZ], ctx)?;
1733 DASADD(HANDLE, (3 * NV), VRTCES.as_slice(), ctx)?;
1734
1735 //
1736 // Next add the integer data to the segment.
1737 //
1738 // NV is the number of vertices.
1739 // NP is the number of plates.
1740 // NVXTOT is the number of voxels in the spatial index.
1741 //
1742 DASADI(HANDLE, 1, &[NV], ctx)?;
1743 DASADI(HANDLE, 1, &[NP], ctx)?;
1744 DASADI(HANDLE, 1, &[NVXTOT], ctx)?;
1745 DASADI(HANDLE, 3, VGREXT.as_slice(), ctx)?;
1746 DASADI(HANDLE, 1, &[CGRSCL], ctx)?;
1747 DASADI(HANDLE, 1, &[VOXNPT], ctx)?;
1748 DASADI(HANDLE, 1, &[VOXNPL], ctx)?;
1749 DASADI(HANDLE, 1, &[VTXNPL], ctx)?;
1750 DASADI(HANDLE, (3 * NP), PLATES.as_slice(), ctx)?;
1751 DASADI(HANDLE, VOXNPT, SPAIXI.subarray(PVOXPT), ctx)?;
1752 DASADI(HANDLE, VOXNPL, SPAIXI.subarray(PVOXPL), ctx)?;
1753 DASADI(HANDLE, NV, SPAIXI.subarray(PVTXPT), ctx)?;
1754 DASADI(HANDLE, VTXNPL, SPAIXI.subarray(PVTXPL), ctx)?;
1755 DASADI(HANDLE, NCGR, SPAIXI.subarray(SICGRD), ctx)?;
1756
1757 //
1758 // End the segment.
1759 //
1760 DLAENS(HANDLE, ctx)?;
1761
1762 CHKOUT(b"DSKW02", ctx)?;
1763 Ok(())
1764}