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}