FreeWRL/FreeX3D  3.0.0
Component_Geospatial.c
1 /*
2 
3 
4 X3D Geospatial Component
5 
6 */
7 
8 
9 /****************************************************************************
10  This file is part of the FreeWRL/FreeX3D Distribution.
11 
12  Copyright 2009 CRC Canada. (http://www.crc.gc.ca)
13 
14  FreeWRL/FreeX3D is free software: you can redistribute it and/or modify
15  it under the terms of the GNU Lesser Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  FreeWRL/FreeX3D is distributed in the hope that it will be useful,
20  but WITHOUT ANY WARRANTY; without even the implied warranty of
21  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  GNU General Public License for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with FreeWRL/FreeX3D. If not, see <http://www.gnu.org/licenses/>.
26 ****************************************************************************/
27 
28 
29 
30 #include <config.h>
31 #include <system.h>
32 #include <display.h>
33 #include <internal.h>
34 
35 #include <libFreeWRL.h>
36 
37 #include "../vrml_parser/Structs.h"
38 #include "../vrml_parser/CRoutes.h"
39 #include "../main/headers.h"
40 
41 #include "../world_script/fieldSet.h"
42 #include "../x3d_parser/Bindable.h"
43 #include "Collision.h"
44 #include "quaternion.h"
45 #include "Viewer.h"
46 #include "../opengl/Frustum.h"
47 #include "../opengl/Material.h"
48 #include "../opengl/OpenGL_Utils.h"
49 #include "../input/EAIHelpers.h" /* for newASCIIString() */
50 
51 #include "Polyrep.h"
52 #include "LinearAlgebra.h"
53 #include "Component_Shape.h" /* for appearance properties */
54 #include "Component_Geospatial.h"
55 #include "Children.h"
56 #include "../scenegraph/RenderFuncs.h"
57 #include "../ui/common.h"
58 
59 /*
60 Coordinate Conversion algorithms were taken from 2 locations after
61 reading and comprehending the references. The code selected was
62 taken and modified, because the original coders "knew their stuff";
63 any problems with the modified code should be sent to John Stewart.
64 
65 ------
66 References:
67 
68 Jean Meeus "Astronomical Algorithms", 2nd Edition, Chapter 11, page 82
69 
70 "World Geodetic System"
71 http://en.wikipedia.org/wiki/WGS84
72 http://en.wikipedia.org/wiki/Geodetic_system#Conversion
73 
74 "Mathworks Aerospace Blockset"
75 http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/index.html?/access/helpdesk/help/toolbox/aeroblks/geocentrictogeodeticlatitude.html&http://www.google.ca/search?hl=en&q=Geodetic+to+Geocentric+conversion+equation&btnG=Google+Search&meta=
76 
77 "TRANSFORMATION OF GEOCENTRIC TO GEODETIC COORDINATES WITHOUT APPROXIMATIONS"
78 http://www.astro.uni.torun.pl/~kb/Papers/ASS/Geod-ASS.htm
79 
80 "Geodetic Coordinate Conversions"
81 http://www.gmat.unsw.edu.au/snap/gps/clynch_pdfs/coordcvt.pdf
82 
83 "TerrestrialCoordinates.c"
84 http://www.lsc-group.phys.uwm.edu/lal/slug/nightly/doxygen/html/TerrestrialCoordinates_8c.html
85 
86 ------
87 Code Conversions:
88 
89 Geodetic to UTM:
90 UTM to Geodetic:
91  Geo::Coordinates::UTM - Perl extension for Latitiude Longitude conversions.
92  Copyright (c) 2000,2002,2004,2007 by Graham Crookham. All rights reserved.
93 
94 Geocentric to Geodetic:
95 Geodetic to Geocentric:
96  Filename: Gdc_To_Gcc_Converter.java
97  Author: Dan Toms, SRI International
98  Package: GeoTransform <http://www.ai.sri.com/geotransform/>
99  Acknowledgements:
100  The algorithms used in the package were created by Ralph Toms and
101  first appeared as part of the SEDRIS Coordinate Transformation API.
102  These were subsequently modified for this package. This package is
103  not part of the SEDRIS project, and the Java code written for this
104  package has not been certified or tested for correctness by NIMA.
105 
106 
107 *********************************************************************/
108 
109 /* defines used to get a SFVec3d into/outof a function that expects a MFVec3d */
110 #define MF_SF_TEMPS struct Multi_Vec3d mIN; struct Multi_Vec3d mOUT; struct Multi_Vec3d gdCoords;
111 #define FREE_MF_SF_TEMPS FREE_IF_NZ(gdCoords.p); FREE_IF_NZ(mOUT.p); FREE_IF_NZ(mIN.p);
112 
113 
114 #define INIT_MF_FROM_SF(myNode, myField) \
115  mIN.n = 1; \
116  mIN.p = MALLOC(struct SFVec3d *, sizeof (struct SFVec3d)); \
117  mIN.p[0].c[0] = myNode-> myField .c[0];\
118  mIN.p[0].c[1] = myNode-> myField .c[1];\
119  mIN.p[0].c[2] = myNode-> myField .c[2];\
120  mOUT.n=0; mOUT.p = NULL; \
121  gdCoords.n=0; gdCoords.p = NULL;
122 
123 #define MF_FIELD_IN_OUT &mIN, &mOUT, &gdCoords
124 #define COPY_MF_TO_SF(myNode, myField) \
125  myNode-> myField .c[0] = mOUT.p[0].c[0]; \
126  myNode-> myField .c[1] = mOUT.p[0].c[1]; \
127  myNode-> myField .c[2] = mOUT.p[0].c[2]; \
128  FREE_IF_NZ(mIN.p); FREE_IF_NZ(mOUT.p);
129 
130 
131 #define MOVE_TO_ORIGIN(me) GeoMove(X3D_GEOORIGIN(me->geoOrigin), &me->__geoSystem, &mIN, &mOUT, &gdCoords);
132 #define COMPILE_GEOSYSTEM(me) compile_geoSystem (me->_nodeType, &me->geoSystem, &me->__geoSystem);
133 
134 #define RADIANS_PER_DEGREE (double)0.0174532925199432957692
135 #define DEGREES_PER_RADIAN (double)57.2957795130823208768
136 
137 #define ENSURE_SPACE(variableInQuestion) \
138  /* enough room for output? */ \
139  if (variableInQuestion ->n < inCoords->n) { \
140  if (variableInQuestion ->p != NULL) { \
141  FREE_IF_NZ(variableInQuestion->p); \
142  } \
143  variableInQuestion ->p = MALLOC(struct SFVec3d *, sizeof (struct SFVec3d) * inCoords->n); \
144  variableInQuestion ->n = inCoords->n; \
145  }
146 
147 /* for UTM, GC, GD conversions */
148 #define ELEVATION_OUT outc->p[i].c[elevation]
149 #define ELEVATION_IN inc->p[i].c[elevation]
150 #define EASTING_IN inc->p[i].c[easting]
151 #define NORTHING_IN inc->p[i].c[northing]
152 #define UTM_SCALE (double)0.9996
153 #define LATITUDE_OUT outc->p[i].c[latitude]
154 #define LONGITUDE_OUT outc->p[i].c[longitude]
155 #define LATITUDE_IN inc->p[i].c[latitude]
156 #define LONGITUDE_IN inc->p[i].c[longitude]
157 
158 #define GC_X_OUT outc->p[i].c[0]
159 #define GC_Y_OUT outc->p[i].c[1]
160 #define GC_Z_OUT outc->p[i].c[2]
161 
162 /* for Gd_Gc conversions */
163 #define GEOSP_AA_A (double)6377563.396
164 #define GEOSP_AA_F (double)299.3249646
165 #define GEOSP_AM_A (double)6377340.189
166 #define GEOSP_AM_F (double)299.3249646
167 #define GEOSP_AN_A (double)6378160
168 #define GEOSP_AN_F (double)298.25
169 #define GEOSP_BN_A (double)6377483.865
170 #define GEOSP_BN_F (double)299.1528128
171 #define GEOSP_BR_A (double)6377397.155
172 #define GEOSP_BR_F (double)299.1528128
173 #define GEOSP_CC_A (double)6378206.4
174 #define GEOSP_CC_F (double)294.9786982
175 #define GEOSP_CD_A (double)6378249.145
176 #define GEOSP_CD_F (double)293.465
177 #define GEOSP_EA_A (double)6377276.345
178 #define GEOSP_EA_F (double)300.8017
179 #define GEOSP_EB_A (double)6377298.556
180 #define GEOSP_EB_F (double)300.8017
181 #define GEOSP_EC_A (double)6377301.243
182 #define GEOSP_EC_F (double)300.8017
183 #define GEOSP_ED_A (double)6377295.664
184 #define GEOSP_ED_F (double)300.8017
185 #define GEOSP_EE_A (double)6377304.063
186 #define GEOSP_EE_F (double)300.8017
187 #define GEOSP_EF_A (double)6377309.613
188 #define GEOSP_EF_F (double)300.8017
189 #define GEOSP_FA_A (double)6378155
190 #define GEOSP_FA_F (double)298.3
191 #define GEOSP_HE_A (double)6378200
192 #define GEOSP_HE_F (double)298.3
193 #define GEOSP_HO_A (double)6378270
194 #define GEOSP_HO_F (double)297
195 #define GEOSP_ID_A (double)6378160
196 #define GEOSP_ID_F (double)298.247
197 #define GEOSP_IN_A (double)6378388
198 #define GEOSP_IN_F (double)297
199 #define GEOSP_KA_A (double)6378245
200 #define GEOSP_KA_F (double)298.3
201 #define GEOSP_RF_A (double)6378137
202 #define GEOSP_RF_F (double)298.257222101
203 #define GEOSP_SA_A (double)6378160
204 #define GEOSP_SA_F (double)298.25
205 #define GEOSP_WD_A (double)6378135
206 #define GEOSP_WD_F (double)298.26
207 #define GEOSP_WE_A (double)6378137
208 #define GEOSP_WE_F (double)298.257223563
209 
210 #define ELLIPSOID(typ) \
211  case typ: Gd_Gc(inCoords,outCoords,typ##_A, typ##_F,geoSystem->p[3], geoSystem->p[4]); break;
212 
213 #define UTM_ELLIPSOID(typ) \
214  case typ: Utm_Gd (inCoords, gdCoords, typ##_A, typ##_F, geoSystem->p[3], geoSystem->p[2], TRUE); \
215  Gd_Gc(gdCoords,outCoords,typ##_A, typ##_F, geoSystem->p[3], geoSystem->p[4]); break;
216 
217 #define GCC_X gcc->c[0]
218 #define GCC_Y gcc->c[1]
219 #define GCC_Z gcc->c[2]
220 #define GDC_LAT gdc->c[0]
221 #define GDC_LON gdc->c[1]
222 #define GDC_ELE gdc->c[2]
223 
224 
225 #define INITIALIZE_GEOSPATIAL(me) \
226  initializeGeospatial((struct X3D_GeoOrigin **) &me->geoOrigin);
227 
228 #define CONVERT_BACK_TO_GD_OR_UTM(thisField) \
229  /* compileGeosystem - encode the return value such that srf->p[x] is... \
230  0: spatial reference frame (GEOSP_UTM, GEOSP_GC, GEOSP_GD); \
231  1: spatial coordinates (defaults to GEOSP_WE) \
232  2: UTM zone number, 1..60. INT_ID_UNDEFINED = not specified \
233  3: UTM: if "S" - value is FALSE, not S, value is TRUE \
234  GD: if "latitude_first" TRUE, if "longitude_first", FALSE \
235  GC: if "northing_first" TRUE, if "easting_first", FALSE */ \
236  \
237  /* do we need to change this from a GCC? */ \
238  if (node->__geoSystem.n != 0) { /* do we have a GeoSystem specified?? if not, dont do this! */ \
239  struct SFVec3d gdCoords; \
240  \
241  if (node->__geoSystem.p[0] != GEOSP_GC) { \
242  /* have to convert to GD or UTM. Go to GD first */ \
243  bool dugsInterpretationOfSpecs = true; \
244  if(dugsInterpretationOfSpecs) \
245  { \
246  retractOrigin((struct X3D_GeoOrigin *)node->geoOrigin, \
247  &thisField); \
248  }else{ \
249  if (Viewer()->GeoSpatialNode != NULL) { \
250  retractOrigin((struct X3D_GeoOrigin *)Viewer()->GeoSpatialNode->geoOrigin, \
251  &thisField); \
252  } \
253  } \
254  \
255  /* printf ("changed retracted, %lf %lf %lf\n", thisField.c[0], thisField.c[1], thisField.c[2]); */ \
256  \
257  /* now, convert to a GDC */ \
258  gccToGdc (&thisField, &gdCoords); \
259  memcpy (&thisField, &gdCoords, sizeof (struct SFVec3d)); \
260  \
261  /* printf ("changed as a GDC, %lf %lf %lf\n", thisField.c[0], thisField.c[1], thisField.c[2]); */ \
262  \
263  /* is this a GD? if so, go no further */ \
264  if (node->__geoSystem.p[0] == GEOSP_GD) { \
265  /* do we need to flip lat and lon? */ \
266  if (!(node->__geoSystem.p[3])) { \
267  double tmp; \
268  tmp = thisField.c[0]; \
269  thisField.c[0] = thisField.c[1]; \
270  thisField.c[1] = tmp; \
271  } \
272  \
273  } else { \
274  /* convert this to UTM */ \
275  int zone; \
276  double easting; \
277  double northing; \
278  \
279  /* get the zone from the geoSystem; if undefined, we will calculate */ \
280  zone = node->__geoSystem.p[2]; \
281  gdToUtm(thisField.c[0], \
282  thisField.c[1], \
283  &zone, &easting, &northing); \
284  \
285  thisField.c[0] = northing; \
286  thisField.c[1] = easting; \
287  \
288  /* printf ("changed as a UTM, %lf %lf %lf\n", thisField[0], thisField[1], thisField[2]); */ \
289  } \
290  } \
291  }
292 
293 //int geoLodLevel = 0;
294 
295 static int gcToGdInit = FALSE;
296 
297 static void compile_geoSystem (int nodeType, struct Multi_String *args, struct Multi_Int32 *srf);
298 static void moveCoords(struct Multi_Int32*, struct Multi_Vec3d *, struct Multi_Vec3d *, struct Multi_Vec3d *);
299 static void Gd_Gc (struct Multi_Vec3d *, struct Multi_Vec3d *, double, double, int, int);
300 static void gccToGdc (struct SFVec3d *, struct SFVec3d *);
301 void calculateViewingSpeed(void);
302 
303 /* for converting from GC to GD */
304 static double A, F, C, A2, C2, Eps2, Eps21, Eps25, C254, C2DA, CEE,
305  CE2, CEEps2, TwoCEE, tem, ARat1, ARat2, BRat1, BRat2, B1,B2,B3,B4,B5;
306 
307 typedef struct pComponent_Geospatial{
308  int geoLodLevel;// = 0;
309 
310 }* ppComponent_Geospatial;
311 void *Component_Geospatial_constructor(){
312  void *v = MALLOCV(sizeof(struct pComponent_Geospatial));
313  memset(v,0,sizeof(struct pComponent_Geospatial));
314  return v;
315 }
316 void Component_Geospatial_init(struct tComponent_Geospatial *t){
317  //public
318  //private
319  t->prv = Component_Geospatial_constructor();
320  {
321  ppComponent_Geospatial p = (ppComponent_Geospatial)t->prv;
322  p->geoLodLevel = 0;
323  }
324 }
325 //ppComponent_Geospatial p = (ppComponent_Geospatial)gglobal()->Component_Geospatial.prv;
326 
327 // http://www.colorado.edu/geography/gcraft/notes/datum/geoid84.html
328 char geoid[][36] = {
329 //-180 longitude ---------------------------------------------------------- + 170 longitude
330 {13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13,13}, //+90 north pole
331 {3,1,-2,-3,-3,-3,-1,3,1,5,9,11,19,27,31,34,33,34,33,34,28,23,17,13,9,4,4,1,-2,-2,0,2,3,2,1,1}, //+80
332 {2,2,1,-1,-3,-7,-14,-24,-27,-25,-19,3,24,37,47,60,61,58,51,43,29,20,12,5,-2,-10,-14,-12,-10,-14,-12,-6,-2,3,6,4}, //+70
333 {2,9,17,10,13,1,-14,-30,-39,-46,-42,-21,6,29,49,65,60,57,47,41,21,18,14,7,-3,-22,-29,-32,-32,-26,-15,-2,13,17,19,6}, //+60
334 {-8,8,8,1,-11,-19,-16,-18,-22,-35,-40,-26,-12,24,45,63,62,59,47,48,42,28,12,-10,-19,-33,-43,-42,-43,-29,-2,17,23,22,6,2}, //+50
335 {-12,-10,-13,-20,-31,-34,-21,-16,-26,-34,-33,-35,-26,2,33,59,52,51,52,48,35,40,33,-9,-28,-39,-48,-59,-50,-28,3,23,37,18,-1,-11}, //+40
336 {-7,-5,-8,-15,-28,-40,-42,-29,-22,-26,-32,-51,-40,-17,17,31,34,44,36,28,29,17,12,-20,-15,-40,-33,-34,-34,-28,7,29,43,20,4,-6}, //+30
337 {5,10,7,-7,-23,-39,-47,-34,-9,-10,-20,-45,-48,-32,-9,17,25,31,31,26,15,6,1,-29,-44,-61,-67,-59,-36,-11,21,39,49,39,22,10}, //+20
338 {13,12,11,2,-11,-28,-38,-29,-10,3,1,-11,-41,-42,-16,3,17,33,22,23,2,-3,-7,-36,-59,-90,-95,-63,-24,12,53,60,58,46,36,26}, //+10
339 {22,16,17,13,1,-12,-23,-20,-14,-3,14,10,-15,-27,-18,3,12,20,18,12,-13,-9,-28,-49,-62,-89,-102,-63,-9,33,58,73,74,63,50,32}, //equator
340 {36,22,11,6,-1,-8,-10,-8,-11,-9,1,32,4,-18,-13,-9,4,14,12,13,-2,-14,-25,-32,-38,-60,-75,-63,-26,0,35,52,68,76,64,52}, //-10
341 {51,27,10,0,-9,-11,-5,-2,-3,-1,9,35,20,-5,-6,-5,0,13,17,23,21,8,-9,-10,-11,-20,-40,-47,-45,-25,5,23,45,58,57,63}, //-20
342 {46,22,5,-2,-8,-13,-10,-7,-4,1,9,32,16,4,-8,4,12,15,22,27,34,29,14,15,15,7,-9,-25,-37,-39,-23,-14,15,33,34,45}, //-30
343 {21,6,1,-7,-12,-12,-12,-10,-7,-1,8,23,15,-2,-6,6,21,24,18,26,31,33,39,41,30,24,13,-2,-20,-32,-33,-27,-14,-2,5,20}, //-40
344 {-15,-18,-18,-16,-17,-15,-10,-10,-8,-2,6,14,13,3,3,10,20,27,25,26,34,39,45,45,38,39,28,13,-1,-15,-22,-22,-18,-15,-14,-10}, //-50
345 {-45,-43,-37,-32,-30,-26,-23,-22,-16,-10,-2,10,20,20,21,24,22,17,16,19,25,30,35,35,33,30,27,10,-2,-14,-23,-30,-33,-29,-35,-43}, //-60
346 {-61,-60,-61,-55,-49,-44,-38,-31,-25,-16,-6,1,4,5,4,2,6,12,16,16,17,21,20,26,26,22,16,10,-1,-16,-29,-36,-46,-55,-54,-59}, //-70
347 {-53,-54,-55,-52,-48,-42,-38,-38,-29,-26,-26,-24,-23,-21,-19,-16,-12,-8,-4,-1,1,4,4,6,5,4,2,-6,-15,-24,-33,-40,-48,-50,-53,-52}, //-80
348 {-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30,-30} //-90 south pole
349 };
350 
351 static double geoidCorrection(double latitudeDeg, double longitudeDeg)
352 {
353  //assumes: GD coords in degrees, -180 to +180 range for longitude, -90 to +90 latitude range
354  //returns: what to add to sea level height to get an ellipsoid height
355  // - does not add it - you do that in the caller.
356  // - http://en.wikipedia.org/wiki/Geoid has a global diagram to check the sign/sense of the correction
357  //benefits: allows you to mix GPS and topographic data in the same scene and have good alignment
358  //usage: put 2 GeoOrigins one for GPS and one for local Topographic in your scene, but set them (initially)
359  // both the same, and set them both without the 'WGS84' (geoid) geosystem option.
360  // Apply option 'WGS84' (geoid) option to your topographic data GeoCoordinate GeoSystem (except GeoOrigin).
361  // Then touch up when viewing them in the same scene by adjusting your local topographic geoOrigin height.
362  int il, ip, il1, ip1;
363  double dl, dp, d00, d01, d10, d11, d;
364  //step 1: find the cell indexes
365  il = (int)(longitudeDeg/10.0) + 18 -1; //longitude cell
366  ip = 18 - ((int)(latitudeDeg/10.0) + 9); //latitude cell
367  il1 = il + 1;
368  ip1 = ip - 1;
369  if(il1 > 35) il1 = 0;
370  if(ip1 < 0) ip1 = 0;
371  //step 2: compute the corners of the cell
372  d00 = (float)geoid[ip][il]; //lower left
373  d01 = (float)geoid[ip1][il]; //upper left
374  d10 = (float)geoid[ip][il1]; //lower right
375  d11 = (float)geoid[ip1][il1]; //upper right
376  //step 3: find the position in the cell
377  dl = longitudeDeg - (-180.0 + (float)(il)*10.0);
378  dp = latitudeDeg - (90.0 - (float)(ip)*10.0);
379  //step 4: find the normalized position within the cell range 0-1
380  dl /= 10.0;
381  dp /= 10.0;
382  //step 5: bilinear interpolate from 4 corners to position in cell
383  d = (1.0 - dl)*(1.0 - dp)*d00 + (dl)*(1.0 - dp)*d10 + (1.0-dl)*(dp)*d01 + (dl)*(dp)*d11;
384  // d is how much higher the geoid is from the ellipsoid.
385  d = -d; // to correct a geoid map to ellpsoid, subtract this amount
386  return (double)d;
387 }
388 /* move ourselves BACK to the from the GeoOrigin */
389 static void retractOrigin(struct X3D_GeoOrigin *myGeoOrigin, struct SFVec3d *gcCoords) {
390  if (myGeoOrigin != NULL) {
391  if(myGeoOrigin->rotateYUp == TRUE)
392  {
393  int i;
394  Quaternion rq;
395  struct SFVec3d temp;
396  vrmlrot_to_quaternion(&rq,myGeoOrigin->__rotyup.c[0], myGeoOrigin->__rotyup.c[1], myGeoOrigin->__rotyup.c[2], myGeoOrigin->__rotyup.c[3]);
397  //quaternion_multi_rotation(outxyz,&rq,inxyz,8);
398  quaternion_rotation((struct point_XYZ *)temp.c, &rq, (const struct point_XYZ *)gcCoords->c);
399  for(i=0;i<3;i++)
400  gcCoords->c[i] = temp.c[i];
401  }
402  gcCoords->c[0] += myGeoOrigin->__movedCoords.c[0];
403  gcCoords->c[1] += myGeoOrigin->__movedCoords.c[1];
404  gcCoords->c[2] += myGeoOrigin->__movedCoords.c[2];
405  }
406 }
407 
408 
409 /* convert GD ellipsiod to GC coordinates */
410 static void Gd_Gc (struct Multi_Vec3d *inc, struct Multi_Vec3d *outc, double radius, double eccentricity, int lat_first, int geoid) {
411  int i;
412  double A = radius;
413  double A2 = radius*radius;
414  double F = (double)(1/eccentricity);
415  double C = A*((double)1.0 - F);
416  double C2 = C*C;
417  double Eps2 = F*((double)2.0 - F);
418  double Eps25 = (double) 0.25 * Eps2;
419 
420  int latitude = 0;
421  int longitude = 1;
422  int elevation = 2;
423 
424  double source_lat;
425  double source_lon;
426  double slat;
427  double slat2;
428  double clat;
429  double Rn;
430  double RnPh;
431 
432  if (!lat_first) {
433  printf ("Gd_Gc, NOT lat first\n");
434  latitude = 1; longitude = 0;
435  }
436 
437  /* enough room for output? */
438  if (outc->n < inc->n) {
439  FREE_IF_NZ(outc->p);
440  outc->p = MALLOC(struct SFVec3d *, sizeof (struct SFVec3d) * inc->n);
441  outc->n = inc->n;
442  }
443  #ifdef VERBOSE
444  printf ("Gd_Gc, have n of %d\n",inc->n);
445  #endif
446 
447  for (i=0; i<inc->n; i++) {
448  #ifdef VERBOSE
449  printf ("Gd_Gc, ining lat %lf long %lf ele %lf ",LATITUDE_IN, LONGITUDE_IN, ELEVATION_IN);
450  #endif
451 
452  source_lat = RADIANS_PER_DEGREE * LATITUDE_IN;
453  source_lon = RADIANS_PER_DEGREE * LONGITUDE_IN;
454 
455  #ifdef VERBOSE
456  printf ("Source Latitude %lf Source Longitude %lf\n",source_lat, source_lon);
457  #endif
458 
459  slat = sin(source_lat);
460  slat2 = slat*slat;
461  clat = cos(source_lat);
462 
463  #ifdef VERBOSE
464  printf ("slat %lf slat2 %lf clat %lf\n",slat, slat2, clat);
465  #endif
466 
467 
468  /* square root approximation for Rn */
469  Rn = A / ( (.25 - Eps25 * slat2 + .9999944354799/4) + (.25-Eps25 * slat2)/(.25 - Eps25 * slat2 + .9999944354799/4));
470 
471  RnPh = Rn + ELEVATION_IN;
472  if(geoid)
473  RnPh += geoidCorrection(LATITUDE_IN,LONGITUDE_IN);
474 
475  #ifdef VERBOSE
476  printf ("Rn %lf RnPh %lf\n",Rn, RnPh);
477  #endif
478 
479  GC_X_OUT = RnPh * clat * cos(source_lon);
480  GC_Y_OUT = RnPh * clat * sin(source_lon);
481  GC_Z_OUT = ((C2 / A2) * Rn + ELEVATION_IN) * slat;
482 
483  #ifdef VERBOSE
484  printf ("Gd_Gc, outing x %lf y %lf z %lf\n", GC_X_OUT, GC_Y_OUT, GC_Z_OUT);
485  #endif
486  }
487 }
488 
489 /* convert UTM to GC coordinates by converting to GD as an intermediary step */
490 static void Utm_Gd (struct Multi_Vec3d *inc, struct Multi_Vec3d *outc, double radius, double flatten, int hemisphere_north, int zone, int northing_first) {
491  int i;
492  int northing = 0; /* for determining which input value is northing */
493  int easting = 1; /* for determining which input value is easting */
494  int elevation = 2; /* elevation is always third value, input AND output */
495  int latitude = 0; /* always return latitude as first value */
496  int longitude = 1; /* always return longtitude as second value */
497 
498  /* create the ERM constants. */
499  double F = 1.0/flatten;
500  double Eccentricity = (F) * (2.0-F);
501 
502  double myEasting;
503  double myphi1rad;
504  double myN1;
505  double myT1;
506  double myC1;
507  double myR1;
508  double myD;
509  double Latitude;
510  double Longitude;
511  double longitudeOrigin;
512  double myeccPrimeSquared;
513  double myNorthing;
514  double northingDRCT1;
515  double eccRoot;
516  double calcConstantTerm1;
517  double calcConstantTerm2;
518  double calcConstantTerm3;
519  double calcConstantTerm4;
520 
521  /* is the values specified with an "easting_first?" */
522  if (!northing_first) { northing = 1; easting = 0; }
523 
524  #ifdef VERBOSE
525  if (!northing_first) printf ("UTM to GD, not northing first, flipping norhting and easting\n");
526  #endif
527 
528  #ifdef VERBOSE
529  if (northing_first) printf ("Utm_Gd: northing first\n"); else printf ("Utm_Gd: NOT northing_first\n");
530  if (!hemisphere_north) printf ("Utm_Gd: NOT hemisphere_north\n"); else printf ("Utm_Gd: hemisphere_north\n");
531  #endif
532 
533 
534  /* enough room for output? */
535  if (outc->n < inc->n) {
536  FREE_IF_NZ(outc->p);
537  outc->p = MALLOC(struct SFVec3d *, sizeof (struct SFVec3d) * inc->n);
538  outc->n = inc->n;
539  }
540 
541  /* constants for all UTM vertices */
542  longitudeOrigin = (zone -1) * 6 - 180 + 3;
543  myeccPrimeSquared = Eccentricity/(((double) 1.0) - Eccentricity);
544  eccRoot = (((double)1.0) - sqrt (((double)1.0) - Eccentricity))/
545  (((double)1.0) + sqrt (((double)1.0) - Eccentricity));
546 
547  calcConstantTerm1 = ((double)1.0) -Eccentricity/
548  ((double)4.0) - ((double)3.0) *Eccentricity*Eccentricity/
549  ((double)64.0) -((double)5.0) *Eccentricity*Eccentricity*Eccentricity/((double)256.0);
550 
551  calcConstantTerm2 = ((double)3.0) * eccRoot/((double)2.0) - ((double)27.0) *eccRoot*eccRoot*eccRoot/((double)32.0);
552  calcConstantTerm3 = ((double)21.0) * eccRoot*eccRoot/ ((double)16.0) - ((double)55.0) *eccRoot*eccRoot*eccRoot*eccRoot/ ((double)32.0);
553  calcConstantTerm4 = ((double)151.0) *eccRoot*eccRoot*eccRoot/ ((double)96.0);
554 
555  #ifdef VERBOSE
556  printf ("zone %d\n",zone);
557  printf ("longitudeOrigin %lf\n",longitudeOrigin);
558  printf ("myeccPrimeSquared %lf\n",myeccPrimeSquared);
559  printf ("eccRoot %lf\n",eccRoot);
560  #endif
561 
562  /* go through each vertex specified */
563  for(i=0;i<inc->n;i++) {
564  /* get the values for THIS UTM vertex */
565  ELEVATION_OUT = ELEVATION_IN;
566  myEasting = EASTING_IN - 500000;
567  if (hemisphere_north) myNorthing = NORTHING_IN;
568  else myNorthing = NORTHING_IN - (double)10000000.0;
569 
570  #ifdef VERBOSE
571  printf ("myEasting %lf\n",myEasting);
572  printf ("myNorthing %lf\n",myNorthing);
573  #endif
574 
575 
576  /* scale the northing */
577  myNorthing= myNorthing / UTM_SCALE;
578 
579  northingDRCT1 = myNorthing /(radius * calcConstantTerm1);
580 
581  myphi1rad = northingDRCT1 +
582  calcConstantTerm2 * sin(((double)2.0) *northingDRCT1)+
583  calcConstantTerm3 * sin(((double)4.0) *northingDRCT1)+
584  calcConstantTerm4 * sin(((double)6.0) *northingDRCT1);
585 
586  myN1 = radius/sqrt(((double)1.0) - Eccentricity * sin(myphi1rad) * sin (myphi1rad));
587  myT1 = tan(myphi1rad) * tan(myphi1rad);
588  myC1 = Eccentricity * cos(myphi1rad) * cos (myphi1rad);
589  myR1 = radius * (((double)1.0) - Eccentricity) / pow(((double)1.0) - Eccentricity * sin(myphi1rad) * sin (myphi1rad), 1.5);
590  myD = myEasting/(myN1*UTM_SCALE);
591 
592  Latitude = myphi1rad-(myN1*tan(myphi1rad)/myR1)*
593  (myD*myD/((double)2.0) -
594  (((double)5.0) + ((double)3.0) *myT1+ ((double)10.0) *myC1-
595  ((double)4.0) *myC1*myC1- ((double)9.0) *myeccPrimeSquared)*
596 
597  myD*myD*myD*myD/((double)24.0) +(((double)61.0) +((double)90.0) *
598  myT1+((double)298.0) *myC1+ ((double)45.0) *myT1*myT1-
599  ((double)252.0) * myeccPrimeSquared- ((double)3.0) *myC1*myC1)*myD*myD*myD*myD*myD*myD/((double)720.0));
600 
601 
602  Longitude = (myD-(((double)1.0)+((double)2.0)*myT1+myC1)*myD*myD*myD/((double)6.0)+(((double)5.0) - ((double)2.0) *myC1+
603  ((double)28.0) *myT1-((double)3.0) *myC1*myC1+
604  ((double)8.0) *myeccPrimeSquared+((double)24.0) *myT1*myT1)*myD*myD*myD*myD*myD/120)/cos(myphi1rad);
605 
606 
607 
608  LATITUDE_OUT = Latitude * DEGREES_PER_RADIAN;
609  LONGITUDE_OUT = longitudeOrigin + Longitude * DEGREES_PER_RADIAN;
610 
611 
612  #ifdef VERBOSE
613  /* printf ("myNorthing scaled %lf\n",myNorthing);
614  printf ("northingDRCT1 %lf\n",northingDRCT1);
615  printf ("myphi1rad %lf\n",myphi1rad);
616  printf ("myN1 %lf\n",myN1);
617  printf ("myT1 %lf\n",myT1);
618  printf ("myC1 %lf\n",myC1);
619  printf ("myR1 %lf\n",myR1);
620  printf ("myD %lf\n",myD);
621  printf ("latitude %lf\n",Latitude);
622  printf ("longitude %lf\n",Longitude);
623  */
624  printf ("utmtogd\tnorthing %lf easting %lf ele %lf\n\tlat %lf long %lf ele %lf\n", NORTHING_IN, EASTING_IN, ELEVATION_IN, LATITUDE_OUT, LONGITUDE_OUT, ELEVATION_IN);
625  #endif
626  }
627 }
628 
629 /* take a set of coords, and a geoSystem, and create a set of moved coords */
630 /* we keep around the GD coords because we need them for rotation calculations */
631 /* parameters:
632  geoSystem: compiled geoSystem integer array pointer
633  inCoords: coordinate structure for input coordinates, ANY coordinate type
634  outCoords: area for GC coordinates. Will MALLOC size if required
635  gdCoords: GD coordinates, used for rotation calculations in later stages. WILL MALLOC THIS */
636 
637 static void moveCoords (struct Multi_Int32* geoSystem, struct Multi_Vec3d *inCoords, struct Multi_Vec3d *outCoords, struct Multi_Vec3d *gdCoords) {
638 
639  int i;
640 
641  /* tmpCoords used for UTM coding */
642  gdCoords->n=0; gdCoords->p=NULL;
643 
644  /* make sure the output has enough space for our converted data */
645  ENSURE_SPACE(outCoords)
646  ENSURE_SPACE(gdCoords)
647 
648  /* GD Geosystem - copy coordinates, and convert them to GC */
649  switch (geoSystem->p[0]) {
650  case GEOSP_GD:
651  /* GD_Gd_Gc_convert (inCoords, outCoords); */
652  switch (geoSystem->p[1]) {
653  ELLIPSOID(GEOSP_AA)
654  ELLIPSOID(GEOSP_AM)
655  ELLIPSOID(GEOSP_AN)
656  ELLIPSOID(GEOSP_BN)
657  ELLIPSOID(GEOSP_BR)
658  ELLIPSOID(GEOSP_CC)
659  ELLIPSOID(GEOSP_CD)
660  ELLIPSOID(GEOSP_EA)
661  ELLIPSOID(GEOSP_EB)
662  ELLIPSOID(GEOSP_EC)
663  ELLIPSOID(GEOSP_ED)
664  ELLIPSOID(GEOSP_EE)
665  ELLIPSOID(GEOSP_EF)
666  ELLIPSOID(GEOSP_FA)
667  ELLIPSOID(GEOSP_HE)
668  ELLIPSOID(GEOSP_HO)
669  ELLIPSOID(GEOSP_ID)
670  ELLIPSOID(GEOSP_IN)
671  ELLIPSOID(GEOSP_KA)
672  ELLIPSOID(GEOSP_RF)
673  ELLIPSOID(GEOSP_SA)
674  ELLIPSOID(GEOSP_WD)
675  ELLIPSOID(GEOSP_WE)
676  default: printf ("unknown Gd_Gc: %s\n", stringGEOSPATIALType(geoSystem->p[1]));
677  }
678 
679  /* now, for the GD coord return values; is this in the correct format for calculating
680  rotations? */
681  gdCoords->n = inCoords->n;
682 
683  /* is the GD value NOT the WGS84 ellipsoid? */
684  if (geoSystem->p[1] != GEOSP_WE) {
685  /*no, convert BACK from the GC to GD, WGS84 level for the gd value returns */
686  for (i=0; i<outCoords->n; i++) {
687  gccToGdc (&outCoords->p[i], &gdCoords->p[i]);
688  }
689  } else {
690  /* just copy the coordinates for the GD temporary return */
691  memcpy (gdCoords->p, inCoords->p, sizeof (struct SFVec3d) * inCoords->n);
692  //Q. should geoid correction be added, so gd are in ellipsoid heights like GPS? (vs sea level heights)
693  if(geoSystem->p[4] == TRUE)
694  for(i=0;i<gdCoords->n;i++)
695  gdCoords->p[i].c[2] += geoidCorrection(gdCoords->p[i].c[1-geoSystem->p[3]],gdCoords->p[i].c[geoSystem->p[3]]);
696  }
697  break;
698  case GEOSP_GC:
699  /* an earth-fixed geocentric coord; no conversion required for gc value returns */
700  for (i=0; i< inCoords->n; i++) {
701  outCoords->p[i].c[0] = inCoords->p[i].c[0];
702  outCoords->p[i].c[1] = inCoords->p[i].c[1];
703  outCoords->p[i].c[2] = inCoords->p[i].c[2];
704 
705  /* convert this coord from GC to GD, WGS84 ellipsoid for gd value returns */
706  gccToGdc (&inCoords->p[i], &gdCoords->p[i]);
707  }
708 
709  break;
710  case GEOSP_UTM:
711  /* GD coords will be returned from the conversion process....*/
712  /* first, convert UTM to GC, then GD, then GD to GC */
713  /* see the compileGeosystem function for geoSystem fields */
714  switch (geoSystem->p[1]) {
715  UTM_ELLIPSOID(GEOSP_AA)
716  UTM_ELLIPSOID(GEOSP_AM)
717  UTM_ELLIPSOID(GEOSP_AN)
718  UTM_ELLIPSOID(GEOSP_BN)
719  UTM_ELLIPSOID(GEOSP_BR)
720  UTM_ELLIPSOID(GEOSP_CC)
721  UTM_ELLIPSOID(GEOSP_CD)
722  UTM_ELLIPSOID(GEOSP_EA)
723  UTM_ELLIPSOID(GEOSP_EB)
724  UTM_ELLIPSOID(GEOSP_EC)
725  UTM_ELLIPSOID(GEOSP_ED)
726  UTM_ELLIPSOID(GEOSP_EE)
727  UTM_ELLIPSOID(GEOSP_EF)
728  UTM_ELLIPSOID(GEOSP_FA)
729  UTM_ELLIPSOID(GEOSP_HE)
730  UTM_ELLIPSOID(GEOSP_HO)
731  UTM_ELLIPSOID(GEOSP_ID)
732  UTM_ELLIPSOID(GEOSP_IN)
733  UTM_ELLIPSOID(GEOSP_KA)
734  UTM_ELLIPSOID(GEOSP_RF)
735  UTM_ELLIPSOID(GEOSP_SA)
736  UTM_ELLIPSOID(GEOSP_WD)
737  UTM_ELLIPSOID(GEOSP_WE)
738  default: printf ("unknown Gd_Gc: %s\n", stringGEOSPATIALType(geoSystem->p[1]));
739  }
740  break;
741  default :
742  printf ("incorrect geoSystem field, %s\n",stringGEOSPATIALType(geoSystem->p[0]));
743  return;
744 
745  }
746 }
747 
748 
749 static void initializeGeospatial (struct X3D_GeoOrigin **nodeptr) {
750  MF_SF_TEMPS
751  struct X3D_GeoOrigin *node = NULL;
752 
753  #ifdef VERBOSE
754  printf ("\ninitializing GeoSpatial code nodeptr %u\n",*nodeptr);
755  #endif
756 
757  if (*nodeptr != NULL) {
758  if (X3D_GEOORIGIN(*nodeptr)->_nodeType != NODE_GeoOrigin) {
759  printf ("expected a GeoOrigin node, but got a node of type %s\n",
760  stringNodeType(X3D_GEOORIGIN(*nodeptr)->_nodeType));
761  *nodeptr = NULL;
762  return;
763  } else {
764  /* printf ("um, just setting geoorign to %u\n",(*nodeptr)); */
765  node = X3D_GEOORIGIN(*nodeptr);
766  }
767 
768  /* printf ("initGeoSpatial ich %d ch %d\n",node->_ichange, node->_change); */
769 
770  if NODE_NEEDS_COMPILING {
771  compile_geoSystem (node->_nodeType, &node->geoSystem, &node->__geoSystem);
772  INIT_MF_FROM_SF(node,geoCoords)
773  moveCoords(&node->__geoSystem, MF_FIELD_IN_OUT);
774  COPY_MF_TO_SF(node, __movedCoords)
775 
776  if(node->rotateYUp == TRUE)
777  {
778  struct SFVec4d orient;
779  int i;
780  Quaternion qz,qx,qr;
781 
782  vrmlrot_to_quaternion (&qz,0.0, 0.0, 1.0, RADIANS_PER_DEGREE*((double)90.0 + gdCoords.p[0].c[1]));
783 
784  #ifdef VERBOSE
785  printf ("GeoOrient qz angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",((double)90.0 + gdCoords->c[1]),
786  RADIANS_PER_DEGREE*((double)90.0 + gdCoords->c[1]),qz.x, qz.y, qz.z,qz.w);
787  #endif
788 
789  vrmlrot_to_quaternion (&qx,1.0, 0.0, 0.0, RADIANS_PER_DEGREE*((double)180.0 - gdCoords.p[0].c[0]));
790 
791  #ifdef VERBOSE
792  printf ("GeoOrient qx angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",
793  ((double)180.0 - gdCoords->c[0]), RADIANS_PER_DEGREE*((double)180.0 - gdCoords->c[0]), qx.x, qx.y, qx.z,qx.w);
794  #endif
795 
796  quaternion_add (&qr, &qx, &qz);
797 
798  #ifdef VERBOSE
799  printf ("GeoOrient qr %lf %lf %lf %lf\n",qr.x, qr.y, qr.z,qr.w);
800  #endif
801 
802  quaternion_to_vrmlrot(&qr, &orient.c[0], &orient.c[1], &orient.c[2], &orient.c[3]);
803  for(i=0;i<4;i++)
804  node->__rotyup.c[i] = orient.c[i];
805  }else{
806  int i;
807  for(i=0;i<4;i++)
808  node->__rotyup.c[i] = 0.0;
809  node->__rotyup.c[1] = 1.0;
810  }
811 
812  #ifdef VERBOSE
813  printf ("initializeGeospatial, __movedCoords %lf %lf %lf, ryup %d, geoSystem %d %d %d %d\n",
814  node->__movedCoords.c[0],
815  node->__movedCoords.c[1],
816  node->__movedCoords.c[2],
817  node->rotateYUp,
818  node->__geoSystem.p[0],
819  node->__geoSystem.p[1],
820  node->__geoSystem.p[2],
821  node->__geoSystem.p[3]);
822  printf ("initializeGeospatial, done\n\n");
823  #endif
824 
825  FREE_MF_SF_TEMPS
826  MARK_NODE_COMPILED
827  }
828  }
829 }
830 
831 /* calculate a translation that moves a Geo node to local space */
832 static void GeoMove(struct X3D_GeoOrigin *geoOrigin, struct Multi_Int32* geoSystem, struct Multi_Vec3d *inCoords, struct Multi_Vec3d *outCoords,
833  struct Multi_Vec3d *gdCoords) {
834  int i;
835  struct X3D_GeoOrigin * myOrigin;
836  Quaternion rq;
837 
838  #ifdef VERBOSE
839  printf ("\nstart of GeoMove... %d coords\n",inCoords->n);
840  #endif
841 
842  /* enough room for output? */
843  if (inCoords->n==0) {return;}
844  if (outCoords->n < inCoords->n) {
845  if (outCoords->n!=0) {
846  FREE_IF_NZ(outCoords->p);
847  }
848  outCoords->p = MALLOC(struct SFVec3d *, sizeof (struct SFVec3d) * inCoords->n);
849  outCoords->n = inCoords->n;
850  }
851 
852  /* set out values to 0.0 for now */
853  for (i=0; i<outCoords->n; i++) {
854  outCoords->p[i].c[0] = (double) 0.0; outCoords->p[i].c[1] = (double) 0.0; outCoords->p[i].c[2] = (double) 0.0;
855  }
856 
857  #ifdef VERBOSE
858  for (i=0; i<outCoords->n; i++) {
859  printf ("start of GeoMove, inCoords %d: %lf %lf %lf\n",i, inCoords->p[i].c[0], inCoords->p[i].c[1], inCoords->p[i].c[2]);
860  }
861  #endif
862 
863 
864 
865  /* check the GeoOrigin attached node */
866  myOrigin = NULL;
867  if (geoOrigin != NULL) {
868  if (X3D_GEOORIGIN(geoOrigin)->_nodeType != NODE_GeoOrigin) {
869  ConsoleMessage ("GeoMove, expected a GeoOrigin, found a %s",stringNodeType(X3D_GEOORIGIN(geoOrigin)->_nodeType));
870  printf ("GeoMove, expected a GeoOrigin, found a %s\n",stringNodeType(X3D_GEOORIGIN(geoOrigin)->_nodeType));
871  return;
872  }
873 
874  myOrigin = geoOrigin; /* local one */
875  }
876  /* printf ("GeoMove, using myOrigin %u, passed in geoOrigin %u with vals %lf %lf %lf\n",myOrigin, myOrigin,
877  myOrigin->geoCoords.c[0], myOrigin->geoCoords.c[1], myOrigin->geoCoords.c[2] ); */
878 
879  moveCoords(geoSystem, inCoords, outCoords, gdCoords);
880 
881 
882  for (i=0; i<outCoords->n; i++) {
883 
884  #ifdef VERBOSE
885  printf ("GeoMove, before subtracting origin %lf %lf %lf\n", outCoords->p[i].c[0], outCoords->p[i].c[1], outCoords->p[i].c[2]);
886  if (myOrigin != NULL) printf (" ... origin %lf %lf %lf\n",myOrigin->__movedCoords.c[0], myOrigin->__movedCoords.c[1], myOrigin->__movedCoords.c[2]);
887  #endif
888 
889  if (myOrigin != NULL) {
890  struct SFVec3d temp;
891 
892  outCoords->p[i].c[0] -= myOrigin->__movedCoords.c[0];
893  outCoords->p[i].c[1] -= myOrigin->__movedCoords.c[1];
894  outCoords->p[i].c[2] -= myOrigin->__movedCoords.c[2];
895  if(myOrigin->rotateYUp == TRUE)
896  {
897  if(i==0)
898  {
899  vrmlrot_to_quaternion(&rq,myOrigin->__rotyup.c[0], myOrigin->__rotyup.c[1], myOrigin->__rotyup.c[2], -myOrigin->__rotyup.c[3]);
900  }
901  //quaternion_multi_rotation(outxyz,&rq,inxyz,8);
902  quaternion_rotation((struct point_XYZ *)temp.c, &rq, (const struct point_XYZ *)outCoords->p[i].c);
903  outCoords->p[i].c[0] = temp.c[0];
904  outCoords->p[i].c[1] = temp.c[1];
905  outCoords->p[i].c[2] = temp.c[2];
906  }
907  }
908 
909  #ifdef VERBOSE
910  printf ("GeoMove, after subtracting origin %lf %lf %lf\n", outCoords->p[i].c[0], outCoords->p[i].c[1], outCoords->p[i].c[2]);
911  #endif
912  }
913 }
914 
915 
916 /* for converting BACK to GD from GC */
917 static void initializeGcToGdParams(void) {
918  A = GEOSP_WE_A;
919  F = GEOSP_WE_F;
920 
921  /* Create the ERM constants. */
922  A2 = A * A;
923  F =1/(F);
924  C =(A) * (1-F);
925  C2 = C * C;
926  Eps2 =(F) * (2.0-F);
927  Eps21 =Eps2 - 1.0;
928  Eps25 =.25 * (Eps2);
929  C254 =54.0 * C2;
930 
931  C2DA = C2 / A;
932  CE2 = A2 - C2;
933  tem = CE2 / C2;
934  CEE = Eps2 * Eps2;
935  TwoCEE =2.0 * CEE;
936  CEEps2 =Eps2 * CE2;
937 
938  /* UPPER BOUNDS ON POINT */
939 
940 
941  ARat1 =pow((A + 50005.0),2);
942  ARat2 =(ARat1) / pow((C+50005.0),2);
943 
944  /* LOWER BOUNDS ON POINT */
945 
946  BRat1 =pow((A-10005.0),2);
947  BRat2 =(BRat1) / pow((C-10005.0),2);
948 
949  /* use WE ellipsoid */
950  B1=0.100225438677758E+01;
951  B2=-0.393246903633930E-04;
952  B3=0.241216653453483E+12;
953  B4=0.133733602228679E+14;
954  B5=0.984537701867943E+00;
955  gcToGdInit = TRUE;
956 }
957 
958 
959 /* convert BACK to a GD coordinate, from GC coordinates using WE ellipsoid */
960 static void gccToGdc (struct SFVec3d *gcc, struct SFVec3d *gdc) {
961  double w2,w,z2,testu,testb,top,top2,rr,q,s12,rnn,s1,zp2,wp,wp2,cf,gee,alpha,cl,arg2,p,xarg,r2,r1,ro,
962  s,roe,arg,v,zo;
963 
964  #ifdef VERBOSE
965  printf ("gccToGdc input %lf %lf %lf\n",GCC_X, GCC_Y, GCC_Z);
966  #endif
967 
968 
969  if (!gcToGdInit) initializeGcToGdParams();
970 
971  w2=GCC_X * GCC_X + GCC_Y * GCC_Y;
972  w=sqrt(w2);
973  z2=GCC_Z * GCC_Z;
974 
975  testu=w2 + ARat2 * z2;
976  testb=w2 + BRat2 * z2;
977 
978  if ((testb > BRat1) && (testu < ARat1))
979  {
980 
981  /*POINT IS BETWEEN-10 KIL AND 50 KIL, SO COMPUTE TANGENT LATITUDE */
982 
983  top= GCC_Z * (B1 + (B2 * w2 + B3) /
984  (B4 + w2 * B5 + z2));
985 
986  top2=top*top;
987 
988  rr=top2+w2;
989 
990  q=sqrt(rr);
991 
992  /* ****************************************************************
993 
994  COMPUTE H IN LINE SQUARE ROOT OF 1-EPS2*SIN*SIN. USE SHORT BINOMIAL
995  EXPANSION PLUS ONE ITERATION OF NEWTON'S METHOD FOR SQUARE ROOTS.
996  */
997 
998  s12=top2/rr;
999 
1000  rnn = A / ( (.25 - Eps25*s12 + .9999944354799/4) + (.25-Eps25*s12)/(.25 - Eps25*s12 + .9999944354799/4));
1001  s1=top/q;
1002 
1003  /******************************************************************/
1004 
1005  /* TEST FOR H NEAR POLE. if SIN(¯)**2 <= SIN(45.)**2 THEN NOT NEAR A POLE.*/
1006 
1007  if (s12 < .50)
1008  GDC_ELE = q-rnn;
1009  else
1010  GDC_ELE = GCC_Z / s1 + (Eps21 * rnn);
1011  GDC_LAT = atan(top / w);
1012  GDC_LON = atan2(GCC_Y,GCC_X);
1013  }
1014  /* POINT ABOVE 50 KILOMETERS OR BELOW -10 KILOMETERS */
1015  else /* Do Exact Solution ************ */
1016  {
1017  wp2=GCC_X * GCC_X + GCC_Y * GCC_Y;
1018  zp2=GCC_Z * GCC_Z;
1019  wp=sqrt(wp2);
1020  cf=C254 * zp2;
1021  gee=wp2 - (Eps21 * zp2) - CEEps2;
1022  alpha=cf / (gee*gee);
1023  cl=CEE * wp2 * alpha / gee;
1024  arg2=cl * (cl + 2.0);
1025  s1=1.0 + cl + sqrt(arg2);
1026  s=pow(s1,(1.0/3.0));
1027  p=alpha / (3.0 * pow(( s + (1.0/s) + 1.0),2));
1028  xarg= 1.0 + (TwoCEE * p);
1029  q=sqrt(xarg);
1030  r2= -p * (2.0 * (1.0 - Eps2) * zp2 / ( q * ( 1.0 + q) ) + wp2);
1031  r1=(1.0 + (1.0 / q));
1032  r2 /=A2;
1033 
1034  /* DUE TO PRECISION ERRORS THE ARGUMENT MAY BECOME NEGATIVE IF SO SET THE ARGUMENT TO ZERO.*/
1035 
1036  if (r1+r2 > 0.0)
1037  ro = A * sqrt( .50 * (r1+r2));
1038  else
1039  ro=0.0;
1040 
1041  ro=ro - p * Eps2 * wp / ( 1.0 + q);
1042  //arg0 = pow(( wp - Eps2 * ro),2) + zp2;
1043  roe = Eps2 * ro;
1044  arg = pow(( wp - roe),2) + zp2;
1045  v=sqrt(arg - Eps2 * zp2);
1046  zo=C2DA * GCC_Z / v;
1047  GDC_ELE = sqrt(arg) * (1.0 - C2DA / v);
1048  top=GCC_Z+ tem*zo;
1049  GDC_LAT = atan( top / wp );
1050  GDC_LON =atan2(GCC_Y,GCC_X);
1051  } /* end of Exact solution */
1052 
1053  GDC_LAT *= DEGREES_PER_RADIAN;
1054  GDC_LON *= DEGREES_PER_RADIAN;
1055 #undef VERBOSE
1056 
1057 }
1058 
1059 /* convert a GDC BACK to a UTM coordinate */
1060 static void gdToUtm(double latitude, double longitude, int *zone, double *easting, double *northing) {
1061 #define DEG2RAD (PI/180.00)
1062 #define GEOSP_WE_INV 0.00669438
1063  double lat_radian;
1064  double long_radian;
1065  double myScale;
1066  int longOrigin;
1067  double longOriginradian;
1068  double eccentprime;
1069  double NNN;
1070  double TTT;
1071  double CCC;
1072  double AAA;
1073  double MMM;
1074 
1075  /* calculate the zone number if it is less than zero. If greater than zero, leave alone! */
1076  if (*zone < 0)
1077  *zone = (int) (((longitude + 180.0)/6.0) + 1);
1078 
1079  lat_radian = latitude * DEG2RAD;
1080  long_radian = longitude * DEG2RAD;
1081  myScale = 0.9996;
1082  longOrigin = (*zone - 1)*6 - 180 + 3;
1083  longOriginradian = longOrigin * DEG2RAD;
1084  eccentprime = GEOSP_WE_INV/(1-GEOSP_WE_INV);
1085 
1086  /*
1087  printf ("lat_radian %lf long_radian %lf myScale %lf longOrigin %d longOriginradian %lf eccentprime %lf\n",
1088  lat_radian, long_radian, myScale, longOrigin, longOriginradian, eccentprime);
1089  */
1090 
1091  NNN = GEOSP_WE_A / sqrt(1-GEOSP_WE_INV * sin(lat_radian)*sin(lat_radian));
1092  TTT = tan(lat_radian) * tan(lat_radian);
1093  CCC = eccentprime * cos(lat_radian)*cos(lat_radian);
1094  AAA = cos(lat_radian) * (long_radian - longOriginradian);
1095  MMM = GEOSP_WE_A
1096  * ( ( 1 - GEOSP_WE_INV/4 - 3 * GEOSP_WE_INV * GEOSP_WE_INV/64
1097  - 5 * GEOSP_WE_INV * GEOSP_WE_INV * GEOSP_WE_INV/256
1098  ) * lat_radian
1099  - ( 3 * GEOSP_WE_INV/8 + 3 * GEOSP_WE_INV * GEOSP_WE_INV/32
1100  + 45 * GEOSP_WE_INV * GEOSP_WE_INV * GEOSP_WE_INV/1024
1101  ) * sin(2 * lat_radian)
1102  + ( 15 * GEOSP_WE_INV * GEOSP_WE_INV/256 +
1103  45 * GEOSP_WE_INV * GEOSP_WE_INV * GEOSP_WE_INV/1024
1104  ) * sin(4 * lat_radian)
1105  - ( 35 * GEOSP_WE_INV * GEOSP_WE_INV * GEOSP_WE_INV/3072
1106  ) * sin(6 * lat_radian)
1107  );
1108 
1109  /* printf ("N %lf T %lf C %lf A %lf M %lf\n",NNN,TTT,CCC,AAA,MMM); */
1110 
1111  *easting = myScale*NNN*(AAA+(1-TTT+CCC)*AAA*AAA*AAA/6
1112  + (5-18*TTT+TTT*TTT+72*CCC-58*eccentprime)*AAA*AAA*AAA*AAA*AAA/120)
1113  + 500000.0;
1114 
1115  *northing= myScale * ( MMM + NNN*tan(lat_radian) *
1116  ( AAA*AAA/2+(5-TTT+9*CCC+4*CCC*CCC)*AAA*AAA*AAA*AAA/24 + (61-58*TTT+TTT*TTT+600*CCC-330*eccentprime) * AAA*AAA*AAA*AAA*AAA*AAA/720));
1117 
1118  /*if (latitude < 0) *northing += 10000000.0;*/
1119 
1120  #ifdef VERBOSE
1121  printf ("gdToUtm: lat %lf long %lf zone %d -> easting %lf northing %lf\n",latitude, longitude, *zone,*easting, *northing);
1122  #endif
1123 }
1124 
1125 /* calculate the rotation needed to apply to this position on the GC coordinate location */
1126 static void GeoOrient (struct X3D_Node *geoOrigin, struct Multi_Int32 *geoSystem, struct SFVec3d *gdCoords, struct SFVec4d *orient) {
1127  Quaternion qx;
1128  Quaternion qz;
1129  Quaternion qr;
1130 
1131  orient->c[0] = 0.0;
1132  orient->c[1] = 1.0;
1133  orient->c[2] = 0.0;
1134  orient->c[3] = 0.0;
1135  /* is this a straight GC geoSystem? If so, we do not do any orientation */
1136  if (geoSystem->n > 0) {
1137  if (geoSystem->p[0] == GEOSP_GC) {
1138  #ifdef VERBOSE
1139  printf ("GeoOrient - simple GC, so no orient\n");
1140  #endif
1141  return;
1142  }
1143  }
1144  if(geoOrigin)
1145  {
1146  if(((struct X3D_GeoOrigin*)geoOrigin)->rotateYUp == TRUE) return;
1147  }
1148 
1149  #ifdef VERBOSE
1150  printf ("GeoOrient - gdCoords->c[0,1] is %f %f\n",gdCoords->c[0],gdCoords->c[1]);
1151  #endif
1152 
1153  /* initialize qx and qz */
1154  vrmlrot_to_quaternion (&qz,0.0, 0.0, 1.0, RADIANS_PER_DEGREE*((double)90.0 + gdCoords->c[1]));
1155 
1156  #ifdef VERBOSE
1157  printf ("GeoOrient qz angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",((double)90.0 + gdCoords->c[1]),
1158  RADIANS_PER_DEGREE*((double)90.0 + gdCoords->c[1]),qz.x, qz.y, qz.z,qz.w);
1159  #endif
1160 
1161  vrmlrot_to_quaternion (&qx,1.0, 0.0, 0.0, RADIANS_PER_DEGREE*((double)180.0 - gdCoords->c[0]));
1162 
1163  #ifdef VERBOSE
1164  printf ("GeoOrient qx angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",
1165  ((double)180.0 - gdCoords->c[0]), RADIANS_PER_DEGREE*((double)180.0 - gdCoords->c[0]), qx.x, qx.y, qx.z,qx.w);
1166  #endif
1167 
1168  quaternion_add (&qr, &qx, &qz);
1169 
1170  #ifdef VERBOSE
1171  printf ("GeoOrient qr %lf %lf %lf %lf\n",qr.x, qr.y, qr.z,qr.w);
1172  #endif
1173 
1174  quaternion_to_vrmlrot(&qr, &orient->c[0], &orient->c[1], &orient->c[2], &orient->c[3]);
1175 
1176  #ifdef VERBOSE
1177  printf ("GeoOrient rotation %lf %lf %lf %lf\n",orient->c[0], orient->c[1], orient->c[2], orient->c[3]);
1178  #endif
1179 }
1180 
1181 
1182 /* compileGeosystem - encode the return value such that srf->p[x] is...
1183  0: spatial reference frame (GEOSP_UTM, GEOSP_GC, GEOSP_GD);
1184  1: spatial coordinates (defaults to GEOSP_WE)
1185  2: UTM zone number, 1..60. INT_ID_UNDEFINED = not specified
1186  3: UTM: if "S" - value is FALSE, not S, value is TRUE
1187  GD: if "latitude_first" TRUE, if "longitude_first", FALSE
1188  GC: if "northing_first" TRUE, if "easting_first", FALSE
1189  4: GD: TRUE if heights are relative to sea level WGS84
1190  (fw converts to ellipsoid heights using geoid correction, allowing
1191  GPS and topographic data to be blended in the same scene)
1192  */
1193 
1194 static void compile_geoSystem (int nodeType, struct Multi_String *args, struct Multi_Int32 *srf) {
1195  int i;
1196  indexT this_srf = INT_ID_UNDEFINED;
1197  indexT this_srf_ind = INT_ID_UNDEFINED;
1198 
1199  #ifdef VERBOSE
1200  printf ("start of compile_geoSystem\n");
1201  #endif
1202 
1203  /* malloc the area required for internal settings, if required */
1204  if (srf->p==NULL) {
1205  srf->n=5;
1206  srf->p=MALLOC(int *, sizeof(int) * 5);
1207  }
1208 
1209  /* set these as defaults */
1210  srf->p[0] = GEOSP_GD;
1211  srf->p[1] = GEOSP_WE;
1212  srf->p[2] = INT_ID_UNDEFINED;
1213  srf->p[3] = TRUE;
1214  srf->p[4] = FALSE; //geoid - not GC, just GD/UTM
1215 
1216  /* if nothing specified, we just use these defaults */
1217  if (args->n==0) return;
1218 
1219  /* first go through, and find the Spatial Reference Frame, GD, UTM, or GC */
1220  for (i=0; i<args->n; i++) {
1221  /* printf ("geoSystem args %d %s\n",i, args->p[i]->strptr); */
1222  indexT tc = findFieldInGEOSPATIAL(args->p[i]->strptr);
1223 
1224  if ((tc == GEOSP_GD) || (tc == GEOSP_GDC)) {
1225  this_srf = GEOSP_GD;
1226  this_srf_ind = i;
1227  } else if ((tc == GEOSP_GC) || (tc == GEOSP_GCC)) {
1228  this_srf = GEOSP_GC;
1229  this_srf_ind = i;
1230  } else if (tc == GEOSP_UTM) {
1231  this_srf = GEOSP_UTM;
1232  this_srf_ind = i;
1233  }
1234  }
1235 
1236  /* did we find a GC, GD, or UTM? */
1237  if (this_srf == INT_ID_UNDEFINED) {
1238  ConsoleMessage ("geoSystem in node %s, must have GC, GD or UTM",stringNodeType(nodeType));
1239  return;
1240  }
1241 
1242  srf->p[0] = (int) this_srf;
1243  /* go through and ensure that we have the correct parameters for this spatial reference frame */
1244  if (this_srf == GEOSP_GC) {
1245  /* possible parameter: GC: if "northing_first" TRUE, if "easting_first", FALSE */
1246  srf->p[1] = INT_ID_UNDEFINED;
1247  for (i=0; i<args->n; i++) {
1248  if (strcmp("northing_first",args->p[i]->strptr) == 0) { srf->p[3] = TRUE;
1249  } else if (strcmp("easting_first",args->p[i]->strptr) == 0) { srf->p[3] = FALSE;
1250  } else if (i!=this_srf_ind) ConsoleMessage ("geoSystem GC parameter %s not allowed geospatial coordinates",args->p[i]->strptr);
1251  }
1252  } else if (this_srf == GEOSP_GD) {
1253  srf->p[1] = GEOSP_WE;
1254  /* possible parameters: ellipsoid, gets put into element 1.
1255  if "latitude_first" TRUE, if "longitude_first", FALSE */
1256 
1257  /* is there an optional argument? */
1258  for (i=0; i<args->n; i++) {
1259  /* printf ("geosp_gd, ind %d i am %d string %s\n",i, this_srf_ind,args->p[i]->strptr); */
1260  if (strcmp("latitude_first", args->p[i]->strptr) == 0) {
1261  srf->p[3] = TRUE;
1262  } else if (strcmp("longitude_first", args->p[i]->strptr) == 0) {
1263  srf->p[3] = FALSE;
1264  } else if(strcmp ("WGS84",args->p[i]->strptr) == 0){
1265  srf->p[4] = TRUE; //geoid
1266  } else {
1267  if (i!= this_srf_ind) {
1268  indexT tc = findFieldInGEOSPATIAL(args->p[i]->strptr);
1269  switch (tc) {
1270  case INT_ID_UNDEFINED:
1271  case GEOSP_GC:
1272  case GEOSP_GCC:
1273  case GEOSP_GD:
1274  case GEOSP_GDC:
1275  case GEOSP_UTM:
1276  ConsoleMessage("expected valid GC parameter in node %s",stringNodeType(nodeType));
1277  srf->p[1] = GEOSP_WE;
1278  break;
1279 
1280  default:
1281  srf->p[1] = (int) tc;
1282  }
1283  }
1284  }
1285  }
1286  } else {
1287  /* this must be UTM */
1288  /* encode the return value such that srf->p[x] is...
1289  0: spatial reference frame (GEOSP_UTM, GEOSP_GC, GEOSP_GD);
1290  1: spatial coordinates (defaults to GEOSP_WE)
1291  2: UTM zone number, 1..60. INT_ID_UNDEFINED = not specified
1292  3: UTM: if "S" - value is FALSE, not S, value is TRUE */
1293  /* first go through, and find the Spatial Reference Frame, GD, UTM, or GC */
1294  for (i=0; i<args->n; i++) {
1295  if (i != this_srf_ind) {
1296  if (strcmp ("S",args->p[i]->strptr) == 0) {
1297  srf->p[3] = FALSE;
1298  } else if (strcmp ("N",args->p[i]->strptr) == 0) {
1299  srf->p[3] = TRUE; // default
1300  } else if (args->p[i]->strptr[0] == 'Z') {
1301  int zone = -1;
1302  sscanf(args->p[i]->strptr,"Z%d",&zone);
1303  /* printf ("zone found as %d\n",zone); */
1304  srf->p[2] = zone;
1305  } else if(strcmp ("WGS84",args->p[i]->strptr) == 0){
1306  srf->p[4] = TRUE; //geoid
1307  } else {
1308  indexT tc = findFieldInGEOSPATIAL(args->p[i]->strptr);
1309  switch (tc) {
1310  case INT_ID_UNDEFINED:
1311  case GEOSP_GC:
1312  case GEOSP_GCC:
1313  case GEOSP_GD:
1314  case GEOSP_GDC:
1315  case GEOSP_UTM:
1316  ConsoleMessage("expected valid UTM Ellipsoid parameter in node %s",stringNodeType(nodeType));
1317  srf->p[1] = GEOSP_WE;
1318  break;
1319 
1320  default:
1321  srf->p[1] = (int)tc;
1322  }
1323  }
1324  }
1325 
1326  }
1327  }
1328  #ifdef VERBOSE
1329  printf ("printf done compileGeoSystem\n");
1330  #endif
1331 
1332 }
1333 
1334 /************************************************************************/
1335 void compile_GeoCoordinate (struct X3D_GeoCoordinate * node) {
1336  MF_SF_TEMPS
1337  int i;
1338 
1339  #ifdef VERBOSE
1340  printf ("compiling GeoCoordinate\n");
1341  #endif
1342 
1343  /* standard MACROS expect specific field names */
1344  mIN = node->point;
1345  mOUT.p = NULL; mOUT.n = 0;
1346 
1347 
1348  INITIALIZE_GEOSPATIAL(node)
1349  COMPILE_GEOSYSTEM(node)
1350  MOVE_TO_ORIGIN(node)
1351 
1352  /* convert the doubles down to floats, because coords are used as floats in FreeWRL. */
1353  FREE_IF_NZ(node->__movedCoords.p);
1354  node->__movedCoords.p = MALLOC (struct SFVec3f *, sizeof (struct SFVec3f) * mOUT.n);
1355  for (i=0; i<mOUT.n; i++) {
1356  node->__movedCoords.p[i].c[0] = (float) mOUT.p[i].c[0];
1357  node->__movedCoords.p[i].c[1] = (float) mOUT.p[i].c[1];
1358  node->__movedCoords.p[i].c[2] = (float) mOUT.p[i].c[2];
1359  #ifdef VERBOSE
1360  printf ("coord %d now is %f %f %f\n", i, node->__movedCoords.p[i].c[0],node->__movedCoords.p[i].c[1],node->__movedCoords.p[i].c[2]);
1361  #endif
1362  }
1363  node->__movedCoords.n = mOUT.n;
1364 
1365  FREE_IF_NZ(gdCoords.p);
1366  FREE_IF_NZ(mOUT.p);
1367  MARK_NODE_COMPILED
1368 
1369  /* events */
1370  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoCoordinate, metadata)) */
1371 }
1372 
1373 
1374 /************************************************************************/
1375 /* GeoElevationGrid */
1376 /************************************************************************/
1377 
1378 /* check validity of ElevationGrid fields */
1379 int checkX3DGeoElevationGridFields (struct X3D_GeoElevationGrid *node, float **points, int *npoints) {
1380  MF_SF_TEMPS
1381  int i,j;
1382  int nx;
1383  double xSp;
1384  int nz;
1385  double zSp;
1386  double *height;
1387  int ntri;
1388  int nh;
1389  struct X3D_PolyRep *rep;
1390  float *newpoints;
1391  int nquads;
1392  int *cindexptr;
1393  float *texcoord = NULL;
1394  double myHeightAboveEllip = 0.0;
1395  int mySRF = 0;
1396 
1397  nx = node->xDimension;
1398  xSp = node->xSpacing;
1399  nz = node->zDimension;
1400  zSp = node->zSpacing;
1401  height = node->height.p;
1402  nh = node->height.n;
1403 
1404  COMPILE_GEOSYSTEM(node)
1405  /* various values for converting to GD/UTM, etc */
1406  if (node->__geoSystem.n != 0) {
1407  mySRF = node->__geoSystem.p[0];
1408  /* NOTE - DO NOT DO THIS CALCULATION - it is added in later
1409  myHeightAboveEllip = getEllipsoidRadius(node->__geoSystem.p[1]);
1410  */
1411  }
1412 
1413  rep = node->_intern;
1414 
1415  /* work out how many triangles/quads we will have */
1416  ntri = (nx && nz ? 2 * (nx-1) * (nz-1) : 0);
1417  nquads = ntri/2;
1418 
1419  /* check validity of input fields */
1420  if(nh != nx * nz) {
1421  if (nh > nx * nz) {
1422  printf ("GeoElevationgrid: warning: x,y vs. height: %d * %d ne %d:\n", nx,nz,nh);
1423  } else {
1424  printf ("GeoElevationgrid: error: x,y vs. height: %d * %d ne %d:\n", nx,nz,nh);
1425  return FALSE;
1426  }
1427  }
1428 
1429  /* do we have any triangles? */
1430  if ((nx < 2) || (nz < 2)) {
1431  printf ("GeoElevationGrid: xDimension and zDimension less than 2 %d %d\n", nx,nz);
1432  return FALSE;
1433  }
1434 
1435  //printf ("checkX3DGeoElevationGrid - node->texCoord %p\n",node->texCoord);
1436 
1437 
1438  /* any texture coordinates passed in? if so, DO NOT generate any texture coords here. */
1439  if (!(node->texCoord)) {
1440  /* allocate memory for texture coords */
1441  FREE_IF_NZ(rep->GeneratedTexCoords[0]);
1442 
1443  /* 6 vertices per quad each vertex has a 2-float tex coord mapping */
1444  texcoord = rep->GeneratedTexCoords[0] = MALLOC (float *, sizeof (float) * nquads * 12);
1445 
1446  rep->tcindex=0; /* we will generate our own mapping */
1447  }
1448 
1449  /* make up points array */
1450  /* a point is a vertex and consists of 3 floats (x,y,z) */
1451  newpoints = MALLOC (float *, sizeof (float) * nz * nx * 3);
1452 
1453  FREE_IF_NZ(rep->actualCoord);
1454  rep->actualCoord = (float *)newpoints;
1455 
1456  /* make up coord index */
1457  if (node->_coordIndex.n > 0) {FREE_IF_NZ(node->_coordIndex.p);}
1458  node->_coordIndex.p = MALLOC (int *, sizeof(int) * nquads * 5);
1459  cindexptr = node->_coordIndex.p;
1460 
1461  node->_coordIndex.n = nquads * 5;
1462  /* return the newpoints array to the caller */
1463  *points = newpoints;
1464  *npoints = node->_coordIndex.n;
1465 
1466  #ifdef VERBOSE
1467  printf ("coordindex:\n");
1468  #endif
1469 
1470  /* ElevationGrids go 1 - 2 - 3 - 4 we go 1 - 4 - 3 - 2 */
1471  //printf ("GeoElevationGrids, nz %d, nx %d\n",nz,nx);
1472 
1473  for (j = 0; j < (nz -1); j++) {
1474  for (i=0; i < (nx-1) ; i++) {
1475  #ifdef VERBOSE
1476  printf (" %d %d %d %d %d\n", j*nx+i, j*nx+i+nx, j*nx+i+nx+1, j*nx+i+1, -1);
1477  #endif
1478 
1479 #ifdef WINDING_ELEVATIONGRID
1480  *cindexptr = j*nx+i; cindexptr++; /* 1 */
1481  *cindexptr = j*nx+i+nx; cindexptr++; /* 2 */
1482  *cindexptr = j*nx+i+nx+1; cindexptr++; /* 3 */
1483  *cindexptr = j*nx+i+1; cindexptr++; /* 4 */
1484  *cindexptr = -1; cindexptr++;
1485 #else
1486  *cindexptr = j*nx+i; cindexptr++; /* 1 */
1487  *cindexptr = j*nx+i+1; cindexptr++; /* 4 */
1488  *cindexptr = j*nx+i+nx+1; cindexptr++; /* 3 */
1489  *cindexptr = j*nx+i+nx; cindexptr++; /* 2 */
1490  *cindexptr = -1; cindexptr++;
1491 #endif
1492 
1493  }
1494  }
1495 
1496  /* tex coords These need to be streamed now; that means for each quad, each vertex needs its tex coords. */
1497  /* if the texCoord node exists, let render_TextureCoordinate (or whatever the node is) do our work for us */
1498  if (!(node->texCoord)) {
1499  //printf ("geoelevationgrid, doing %d x %d texture coords; tcoord %p\n",nz-1,nx-1,texcoord);
1500 
1501  for (j = 0; j < (nz -1); j++) {
1502  for (i=0; i < (nx-1) ; i++) {
1503  /* first triangle, 3 vertexes */
1504 #ifdef WINDING_ELEVATIONGRID
1505  /* first tri */
1506 /* 1 */ *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
1507  *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
1508 
1509 /* 2 */ *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
1510  *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
1511 
1512 /* 3 */ *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
1513  *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
1514 
1515  /* second tri */
1516 /* 1 */ *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
1517  *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
1518 
1519 /* 3 */ *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
1520  *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
1521 
1522 /* 4 */ *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
1523  *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
1524 #else
1525  /* first tri */
1526 /* 1 */ *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
1527  *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
1528 
1529 /* 4 */ *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
1530  *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
1531 
1532 /* 3 */ *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
1533  *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
1534 
1535  /* second tri */
1536 /* 1 */ *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
1537  *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
1538 
1539 /* 3 */ *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
1540  *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
1541 
1542 /* 2 */ *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
1543  *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
1544 
1545 #endif
1546  }
1547  }
1548  //for (i=0; i<10; i++) printf ("geoele tc %d is %f\n",i,rep->GeneratedTexCoords[i]);
1549  }
1550 
1551  /* Render_Polyrep will use this number of triangles */
1552  rep->ntri = ntri;
1553 
1554  /* initialize arrays used for passing values into/out of the MOVE_TO_ORIGIN(node) values */
1555  mIN.n = nx * nz;
1556  mIN.p = MALLOC (struct SFVec3d *, sizeof (struct SFVec3d) * mIN.n);
1557 
1558  mOUT.n=0; mOUT.p = NULL;
1559  gdCoords.n=0; gdCoords.p = NULL;
1560 
1561  /* make up a series of points, then go and convert them to local coords */
1562  for (j=0; j<nz; j++) {
1563  for (i=0; i < nx; i++) {
1564 
1565  #ifdef VERBOSE
1566  printf (" %lf %lf %lf # (hei ind %d) point [%d, %d]\n",
1567  xSp * i,
1568  height[i+(j*nx)] * ((double)node->yScale),
1569  zSp * j,
1570  i+(j*nx), i,j);
1571  #endif
1572 
1573 
1574  /* Make up a new vertex. Add the geoGridOrigin to every point */
1575 
1576  if ((mySRF == GEOSP_GD) || (mySRF == GEOSP_UTM)) {
1577  /* GD - give it to em in Latitude/Longitude/Elevation order */
1578  /* UTM- or give it to em in Northing/Easting/Elevation order */
1579  /* latitude - range of -90 to +90 */
1580  mIN.p[i+(j*nx)].c[0] = zSp * j + node->geoGridOrigin.c[0];
1581 
1582  /* longitude - range -180 to +180, or 0 to 360 */
1583  mIN.p[i+(j*nx)].c[1] =xSp * i + node->geoGridOrigin.c[1];
1584 
1585  /* elevation, above geoid */
1586  mIN.p[i+(j*nx)].c[2] = (height[i+(j*nx)] *(node->yScale)) + node->geoGridOrigin.c[2]
1587  + myHeightAboveEllip;
1588  } else {
1589  /* nothing quite specified here - what do we really do??? */
1590  mIN.p[i+(j*nx)].c[0] = zSp * j + node->geoGridOrigin.c[0];
1591 
1592  mIN.p[i+(j*nx)].c[1] =xSp * i + node->geoGridOrigin.c[1];
1593 
1594  mIN.p[i+(j*nx)].c[2] = (height[i+(j*nx)] *(node->yScale)) + node->geoGridOrigin.c[2]
1595  + myHeightAboveEllip;
1596 
1597  }
1598  /* printf ("height made up of %lf, geoGridOrigin %lf, myHeightAboveEllip %lf\n",(height[i+(j*nx)] *(node->yScale)),node->geoGridOrigin.c[2], myHeightAboveEllip); */
1599  }
1600  }
1601  #ifdef VERBOSE
1602  printf ("points before moving origin:\n");
1603  for (j=0; j<nz; j++) {
1604  for (i=0; i < nx; i++) {
1605  printf (" %lf %lf %lf # lat/long/height before MOVE, index %d\n",mIN.p[i+(j*nx)].c[0],
1606  mIN.p[i+(j*nx)].c[1],mIN.p[i+(j*nx)].c[2],i+(j*nx));
1607 
1608  }
1609  }
1610  #endif
1611 
1612  /* convert this point to a local coordinate */
1613  MOVE_TO_ORIGIN(node)
1614 
1615  /* copy the resulting array back to the ElevationGrid */
1616 
1617  #ifdef VERBOSE
1618  printf ("points:\n");
1619  #endif
1620 
1621  for (j=0; j<nz; j++) {
1622  for (i=0; i < nx; i++) {
1623  /* copy this coordinate into our ElevationGrid array */
1624  newpoints[0] = (float) mOUT.p[i+(j*nx)].c[0];
1625  newpoints[1] = (float) mOUT.p[i+(j*nx)].c[1];
1626  newpoints[2] = (float) mOUT.p[i+(j*nx)].c[2];
1627 
1628  #ifdef VERBOSE
1629  printf (" %f %f %f # converted, index %d\n",newpoints[0],newpoints[1],newpoints[2],i+(j*nx));
1630  #endif
1631 
1632  newpoints += 3;
1633  }
1634  }
1635  FREE_MF_SF_TEMPS;
1636  return TRUE;
1637 }
1638 
1639 
1640 /* a GeoElevationGrid creates a "real" elevationGrid node as a child for rendering. */
1641 void compile_GeoElevationGrid (struct X3D_GeoElevationGrid * node) {
1642 
1643  #ifdef VERBOSE
1644  printf ("compiling GeoElevationGrid\n");
1645  #endif
1646  printf ("compiling GeoElevationGrid\n");
1647 
1648  INITIALIZE_GEOSPATIAL(node)
1649  COMPILE_GEOSYSTEM(node)
1650  MARK_NODE_COMPILED
1651 
1652  /* events */
1653  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoElevationGrid, metadata)) */
1654 
1655 }
1656 
1657 
1658 void render_GeoElevationGrid (struct X3D_GeoElevationGrid *node) {
1659  INITIALIZE_GEOSPATIAL(node)
1660  COMPILE_POLY_IF_REQUIRED (NULL, NULL, node->color, node->normal, node->texCoord)
1661  CULL_FACE(node->solid)
1662  render_polyrep(node);
1663 }
1664 
1665 /************************************************************************/
1666 /* GeoLocation */
1667 /************************************************************************/
1668 
1669 void compile_GeoLocation (struct X3D_GeoLocation * node) {
1670  // JAS int i;
1671  MF_SF_TEMPS
1672 
1673  #ifdef VERBOSE
1674  printf ("compiling GeoLocation\n");
1675  #endif
1676 
1677  /* work out the position */
1678  INITIALIZE_GEOSPATIAL(node)
1679  COMPILE_GEOSYSTEM(node)
1680  INIT_MF_FROM_SF(node, geoCoords)
1681  MOVE_TO_ORIGIN(node)
1682  COPY_MF_TO_SF(node, __movedCoords)
1683 
1684  /* work out the local orientation */
1685  GeoOrient(node->geoOrigin, &node->__geoSystem, &gdCoords.p[0], &node->__localOrient);
1686 
1687  #ifdef VERBOSE
1688  printf ("compile_GeoLocation, orig coords %lf %lf %lf, moved %lf %lf %lf\n", node->geoCoords.c[0], node->geoCoords.c[1], node->geoCoords.c[2], node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
1689  printf (" rotation is %lf %lf %lf %lf\n",
1690  node->__localOrient.c[0],
1691  node->__localOrient.c[1],
1692  node->__localOrient.c[2],
1693  node->__localOrient.c[3]);
1694  #endif
1695 
1696  /* did the geoCoords change?? */
1697  MARK_SFVEC3D_INOUT_EVENT(node->geoCoords, node->__oldgeoCoords, offsetof (struct X3D_GeoLocation, geoCoords))
1698 
1699  /* how about the children field ?? */
1700  MARK_MFNODE_INOUT_EVENT(node->children, node->__oldChildren, offsetof (struct X3D_GeoLocation, children))
1701 
1702  REINITIALIZE_SORTED_NODES_FIELD(node->children,node->_sortedChildren);
1703  MARK_NODE_COMPILED
1704  FREE_MF_SF_TEMPS
1705 
1706  /* events */
1707  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoLocation, metadata)) */
1708 
1709  INITIALIZE_EXTENT;
1710 
1711  #ifdef VERBOSE
1712  printf ("compiled GeoLocation\n\n");
1713  #endif
1714 }
1715 
1716 void child_GeoLocation (struct X3D_GeoLocation *node) {
1717  CHILDREN_COUNT
1718  //LOCAL_LIGHT_SAVE
1719  INITIALIZE_GEOSPATIAL(node)
1720  COMPILE_IF_REQUIRED
1721 
1722  OCCLUSIONTEST
1723 
1724 
1725  /* {
1726  int x;
1727  struct X3D_Node *xx;
1728 
1729  printf ("child_GeoLocation, this %d \n",node);
1730  for (x=0; x<nc; x++) {
1731  xx = X3D_NODE(node->children.p[x]);
1732  printf (" ch %d type %s dist %f\n",node->children.p[x],stringNodeType(xx->_nodeType),xx->_dist);
1733  }
1734  } */
1735 
1736  /* Check to see if we have to check for collisions for this transform. */
1737 
1738  RETURN_FROM_CHILD_IF_NOT_FOR_ME
1739 
1740  /* do we have a local for a child? */
1741  //LOCAL_LIGHT_CHILDREN(node->children);
1742  prep_sibAffectors((struct X3D_Node*)node,&node->__sibAffectors);
1743 
1744  /* now, just render the non-directionalLight children */
1745 
1746  /* printf ("GeoLocation %d, flags %d, render_sensitive %d\n",
1747  node,node->_renderFlags,render_sensitive); */
1748 
1749  #ifdef CHILDVERBOSE
1750  printf ("GeoLocation - doing normalChildren\n");
1751  #endif
1752 
1753  normalChildren(node->children);
1754 
1755  #ifdef CHILDVERBOSE
1756  printf ("GeoLocation - done normalChildren\n");
1757  #endif
1758 
1759  //LOCAL_LIGHT_OFF
1760  prep_sibAffectors((struct X3D_Node*)node,&node->__sibAffectors);
1761 
1762 }
1763 
1764 /* do transforms, calculate the distance */
1765 void prep_GeoLocation (struct X3D_GeoLocation *node) {
1766  INITIALIZE_GEOSPATIAL(node)
1767  COMPILE_IF_REQUIRED
1768 
1769  /* rendering the viewpoint means doing the inverse transformations in reverse order (while poping stack),
1770  * so we do nothing here in that case -ncoder */
1771 
1772  /* printf ("prep_GeoLocation, render_hier vp %d geom %d light %d sens %d blend %d prox %d col %d\n",
1773  render_vp,render_geom,render_light,render_sensitive,render_blend,render_proximity,render_collision); */
1774 
1775  /* do we have any geometry visible, and are we doing anything with geometry? */
1776  OCCLUSIONTEST
1777 
1778  if(!renderstate()->render_vp) {
1779  FW_GL_PUSH_MATRIX();
1780 
1781  /* TRANSLATION */
1782  FW_GL_TRANSLATE_D(node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
1783 
1784  //printf ("prep_GeoLoc trans to %lf %lf %lf\n",node->__movedCoords.c[0],node->__movedCoords.c[1],node->__movedCoords.c[2]);
1785 
1786  FW_GL_ROTATE_RADIANS(node->__localOrient.c[3], node->__localOrient.c[0],node->__localOrient.c[1],node->__localOrient.c[2]);
1787 
1788  /*
1789  printf ("geoLocation trans %7.4f %7.4f %7.4f\n",node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
1790  printf ("geoLocation rotat %7.4f %7.4f %7.4f %7.4f\n",my_rotation, node->__localOrient.c[0],node->__localOrient.c[1],node->__localOrient.c[2]);
1791  */
1792 
1793  /* did either we or the Viewpoint move since last time? */
1794  RECORD_DISTANCE
1795  }
1796 }
1797 void fin_GeoLocation (struct X3D_GeoLocation *node) {
1798  INITIALIZE_GEOSPATIAL(node)
1799  COMPILE_IF_REQUIRED
1800  OCCLUSIONTEST
1801 
1802  if(!renderstate()->render_vp) {
1803  FW_GL_POP_MATRIX();
1804  } else {
1805  if ((node->_renderFlags & VF_Viewpoint) == VF_Viewpoint) {
1806  FW_GL_ROTATE_RADIANS(-node->__localOrient.c[3], node->__localOrient.c[0],node->__localOrient.c[1],node->__localOrient.c[2]);
1807  FW_GL_TRANSLATE_D(-node->__movedCoords.c[0], -node->__movedCoords.c[1], -node->__movedCoords.c[2]);
1808 
1809  }
1810  }
1811 }
1812 
1813 /************************************************************************/
1814 /* GeoLOD */
1815 /************************************************************************/
1816 void add_node_to_broto_context(struct X3D_Proto *currentContext,struct X3D_Node *node);
1817 
1818 void deleteMallocedFieldValue(int type,union anyVrml *fieldPtr);
1819 void LOAD_CHILD(struct X3D_GeoLOD *node, struct X3D_Node **childNode, struct Multi_String *childUrl) {
1820  /* printf ("start of LOAD_CHILD, url has %d strings\n",node->childUrl.n); */
1821  int i;
1822  if (childUrl->n > 0) {
1823  /* create new inline node, link it in */
1824  if (*childNode == NULL) {
1825  *childNode = createNewX3DNode(NODE_Inline);
1826  if(node->_executionContext)
1827  add_node_to_broto_context(X3D_PROTO(node->_executionContext),X3D_NODE(*childNode));
1828  ADD_PARENT(X3D_NODE(*childNode), X3D_NODE(node));
1829  }
1830  /* copy over the URL from parent */
1831  deleteMallocedFieldValue(FIELDTYPE_MFString,(union anyVrml*)&X3D_INLINE(*childNode)->url);
1832  X3D_INLINE(*childNode)->url.p = MALLOC(struct Uni_String **, sizeof(struct Uni_String)*childUrl->n);
1833  for (i=0; i<childUrl->n; i++) {
1834  /* printf ("copying over url %s\n",node->childUrl.p[i]->strptr); */
1835  X3D_INLINE(*childNode)->url.p[i] = newASCIIString(childUrl->p[i]->strptr);
1836  }
1837  /* printf ("loading, and urlCount is %d\n",node->childUrl.n); */
1838  X3D_INLINE(*childNode)->url.n = childUrl->n;
1839  X3D_INLINE(*childNode)->load = TRUE;
1840  }
1841 }
1842 
1843 #define UNLOAD_CHILD(childNode) \
1844  if (node->childNode != NULL) \
1845  X3D_INLINE(node->childNode)->load = FALSE;
1846 
1847 
1848 static void GeoLODchildren (struct X3D_GeoLOD *node) {
1849  int load = node->__inRange;
1850 
1851  /* lets see if we still have to load this one... */
1852  if (((node->__childloadstatus)==0) && (load)) {
1853  #ifdef VERBOSE
1854  ppComponent_Geospatial p = (ppComponent_Geospatial)gglobal()->Component_Geospatial.prv;
1855 
1856  printf ("GeoLODchildren - have to LOAD_CHILD for node %u (level %d)\n",node,p->geoLodLevel);
1857  #endif
1858 
1859  LOAD_CHILD(node,&node->__child1Node,&node->child1Url);
1860  LOAD_CHILD(node,&node->__child2Node,&node->child2Url);
1861  LOAD_CHILD(node,&node->__child3Node,&node->child3Url);
1862  LOAD_CHILD(node,&node->__child4Node,&node->child4Url);
1863 
1864  //LOAD_CHILD(__child1Node,child1Url)
1865  //LOAD_CHILD(__child2Node,child2Url)
1866  //LOAD_CHILD(__child3Node,child3Url)
1867  //LOAD_CHILD(__child4Node,child4Url)
1868  node->__childloadstatus = 1;
1869  }
1870 }
1871 //void GeoLODchildren1 (struct X3D_GeoLOD *node) {
1872 // GeoLODchildren(node);
1873 //}
1874 static void GeoUnLODchildren (struct X3D_GeoLOD *node) {
1875  int load = node->__inRange;
1876 
1877  if (!(load) && ((node->__childloadstatus) != 0)) {
1878  #ifdef VERBOSE
1879  ppComponent_Geospatial p = (ppComponent_Geospatial)gglobal()->Component_Geospatial.prv;
1880  printf ("GeoLODloadChildren, removing children from node %u level %d\n",node,p->geoLodLevel);
1881  #endif
1882  UNLOAD_CHILD(__child1Node)
1883  UNLOAD_CHILD(__child2Node)
1884  UNLOAD_CHILD(__child3Node)
1885  UNLOAD_CHILD(__child4Node)
1886 
1887  node->__childloadstatus = 0;
1888  }
1889 }
1890 
1891 
1892 static void GeoLODrootUrl (struct X3D_GeoLOD *node) {
1893  int load = node->__inRange == 0; //dug9 it's when you are out of range that you should get the rootnode
1894 
1895  /* lets see if we still have to load this one... */
1896  if (((node->__rooturlloadstatus)==0) && (load)) {
1897  #ifdef VERBOSE
1898  printf ("GeoLODrootUrl - have to LOAD_CHILD for node %u\n",node);
1899  #endif
1900 
1901  LOAD_CHILD(node,&node->__rootUrl, &node->rootUrl);
1902  //LOAD_CHILD(__rootUrl, rootUrl)
1903 
1904  node->__rooturlloadstatus = 1;
1905  }
1906 }
1907 
1908 
1909 static void GeoUnLODrootUrl (struct X3D_GeoLOD *node) {
1910  int load = node->__inRange;
1911 
1912  if (!(load) && ((node->__rooturlloadstatus) != 0)) {
1913  #ifdef VERBOSE
1914  printf ("GeoLODloadChildren, removing rootUrl\n");
1915  #endif
1916  node->__childloadstatus = 0;
1917  }
1918 }
1919 
1920 
1921 
1922 void compile_GeoLOD (struct X3D_GeoLOD * node) {
1923  MF_SF_TEMPS
1924 
1925  #ifdef VERBOSE
1926  printf ("compiling GeoLOD %u\n",node);
1927  #endif
1928 
1929  /* work out the position */
1930  INITIALIZE_GEOSPATIAL(node)
1931  COMPILE_GEOSYSTEM(node)
1932  INIT_MF_FROM_SF(node, center)
1933  MOVE_TO_ORIGIN(node)
1934  COPY_MF_TO_SF(node, __movedCoords)
1935 
1936  #ifdef VERBOSE
1937  printf ("compile_GeoLOD %u, orig coords %lf %lf %lf, moved %lf %lf %lf\n", node, node->center.c[0], node->center.c[1], node->center.c[2], node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
1938 
1939  printf ("children.n %d childurl 1: %u 2: %u 3: %u 4: %u rootUrl: %u rootNode: %d\n",
1940  node->children,
1941  node->child1Url,
1942  node->child2Url,
1943  node->child3Url,
1944  node->child4Url,
1945  node->rootUrl,
1946  node->rootNode.n);
1947  #endif
1948 
1949  MARK_NODE_COMPILED
1950  FREE_MF_SF_TEMPS
1951 
1952  /* events */
1953  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoLOD, metadata)) */
1954 
1955 
1956  #ifdef VERBOSE
1957  printf ("compiled GeoLOD\n\n");
1958  #endif
1959 }
1960 #undef VERBOSE
1961 
1962 
1963 void child_GeoLOD (struct X3D_GeoLOD *node) {
1964  int i;
1965  ppComponent_Geospatial p = (ppComponent_Geospatial)gglobal()->Component_Geospatial.prv;
1966 
1967  INITIALIZE_GEOSPATIAL(node)
1968  COMPILE_IF_REQUIRED
1969 
1970  #ifdef VERBOSE
1971  printf ("child_GeoLOD %u (level %d), renderFlags %x render_hier vp %d geom %d light %d sens %d blend %d prox %d col %d\n",
1972  node,
1973  p->geoLodLevel,
1974  node->_renderFlags,
1975  render_vp,render_geom,render_light,render_sensitive,render_blend,render_proximity,render_collision);
1976  #endif
1977  //ConsoleMessage("glod kids=%d\r",node->children.n);
1978  /* for debugging purposes... */
1979  if (node->__level == -1) node->__level = p->geoLodLevel;
1980  else if (node->__level != p->geoLodLevel) {
1981  printf ("hmmm - GeoLOD %p was level %d, now %d\n",node,node->__level, p->geoLodLevel);
1982  }
1983 
1984  #ifdef VERBOSE
1985  if ( node->__inRange) {
1986  printf ("GeoLOD %u (level %d) closer\n",node,p->geoLodLevel);
1987  } else {
1988  printf ("GeoLOD %u (level %d) farther away\n",node,p->geoLodLevel);
1989  }
1990  #endif
1991 
1992  /* if we are out of range, use the rootNode or rootUrl field */
1993  /* else, use the child1Url through the child4Url fields */
1994  if (!(node->__inRange)) {
1995  /* printf ("GeoLOD, node %u, doing rootNode, rootNode.n = %d\n",node,node->rootNode.n); */
1996  /* do we need to unload children that are no longer needed? */
1997  GeoUnLODchildren (node);
1998 
1999  if (node->rootNode.n != 0) {
2000  for (i=0; i<node->rootNode.n; i++) {
2001  #ifdef VERBOSE
2002  printf ("GeoLOD %u is rendering rootNode %u",node,node->rootNode.p[i]);
2003  if (node->rootNode.p[i]!=NULL) printf (" (%s) ",stringNodeType(X3D_NODE(node->rootNode.p[i])->_nodeType));
2004  printf("\n");
2005  #endif
2006 
2007  render_node (node->rootNode.p[i]);
2008  }
2009  } else if (node->rootUrl.n != 0) {
2010 
2011  /* try and load the root from the rootUrl */
2012  GeoLODrootUrl (node);
2013 
2014  /* render this rootUrl */
2015  if (node->__rootUrl != NULL) {
2016  #ifdef VERBOSE
2017  printf ("GeoLOD %u is rendering rootUrl %u",node,node->__rootUrl);
2018  if (node->__rootUrl != NULL) printf (" (%s) ", stringNodeType(X3D_NODE(node->__rootUrl)->_nodeType));
2019  printf ("\n");
2020  #endif
2021 
2022  render_node (node->__rootUrl);
2023  }
2024 
2025 
2026  }
2027  } else {
2028  p->geoLodLevel++;
2029 
2030  /* go through 4 kids */
2031  GeoLODchildren (node);
2032 
2033  /* get rid of the rootUrl node, if it is loaded */
2034  GeoUnLODrootUrl (node);
2035 
2036  #ifdef VERBOSE
2037  printf ("rendering children at %d, they are: ",p->geoLodLevel);
2038  if (node->child1Url.n>0) printf (" :%s: ",node->child1Url.p[0]->strptr);
2039  if (node->child2Url.n>0) printf (" :%s: ",node->child2Url.p[0]->strptr);
2040  if (node->child3Url.n>0) printf (" :%s: ",node->child3Url.p[0]->strptr);
2041  if (node->child4Url.n>0) printf (" :%s: ",node->child4Url.p[0]->strptr);
2042  printf ("\n");
2043  #endif
2044 
2045  /* render these children */
2046  #ifdef VERBOSE
2047  printf ("GeoLOD %u is rendering children %u ", node, node->__child1Node);
2048  if (node->__child1Node != NULL) printf (" (%s) ",stringNodeType(X3D_NODE(node->__child1Node)->_nodeType));
2049  printf (" %u ", node->__child2Node);
2050  if (node->__child2Node != NULL) printf (" (%s) ",stringNodeType(X3D_NODE(node->__child2Node)->_nodeType));
2051  printf (" %u ", node->__child3Node);
2052  if (node->__child3Node != NULL) printf (" (%s) ",stringNodeType(X3D_NODE(node->__child3Node)->_nodeType));
2053  printf (" %u ", node->__child4Node);
2054  if (node->__child4Node != NULL) printf (" (%s) ",stringNodeType(X3D_NODE(node->__child4Node)->_nodeType));
2055  printf ("\n");
2056  #endif
2057 
2058  if (node->__child1Node != NULL) render_node (node->__child1Node);
2059  if (node->__child2Node != NULL) render_node (node->__child2Node);
2060  if (node->__child3Node != NULL) render_node (node->__child3Node);
2061  if (node->__child4Node != NULL) render_node (node->__child4Node);
2062  p->geoLodLevel--;
2063 
2064  }
2065 }
2066 
2067 /************************************************************************/
2068 /* GeoMetaData */
2069 /************************************************************************/
2070 
2071 void compile_GeoMetadata (struct X3D_GeoMetadata * node) {
2072  #ifdef VERBOSE
2073  printf ("compiling GeoMetadata\n");
2074 
2075  #endif
2076 
2077  MARK_NODE_COMPILED
2078 }
2079 
2080 /************************************************************************/
2081 /* GeoOrigin */
2082 /************************************************************************/
2083 
2084 void compile_GeoOrigin (struct X3D_GeoOrigin * node) {
2085  #ifdef VERBOSE
2086  printf ("compiling GeoOrigin\n");
2087  #endif
2088 
2089  ConsoleMessage ("compiling GeoOrigin\n"); //this doesn't get called - see line 654 in initializeGeospatial()
2090  /* INITIALIZE_GEOSPATIAL */
2091  COMPILE_GEOSYSTEM(node)
2092  {
2093  int i;
2094  for(i=0;i<4;i++)
2095  node->__rotyup.c[0] = 0.0;
2096  node->__rotyup.c[1] = 1.0;
2097  }
2098  MARK_NODE_COMPILED
2099 
2100  /* events */
2101  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoOrigin, metadata)) */
2102  MARK_SFVEC3D_INOUT_EVENT(node->geoCoords,node->__oldgeoCoords,offsetof (struct X3D_GeoOrigin, geoCoords))
2103  //dug9 may 2015 commented out __old.. see also geoViewpoint
2104  //MARK_MFSTRING_INOUT_EVENT(node->geoSystem,node->__oldMFString,offsetof (struct X3D_GeoOrigin, geoSystem))
2105 }
2106 
2107 /************************************************************************/
2108 /* GeoPositionInterpolator */
2109 /************************************************************************/
2110 
2111 void compile_GeoPositionInterpolator (struct X3D_GeoPositionInterpolator * node) {
2112  MF_SF_TEMPS
2113 
2114  #ifdef VERBOSE
2115  printf ("compiling GeoPositionInterpolator\n");
2116  #endif
2117 
2118  /* standard MACROS expect specific field names */
2119  mIN = node->keyValue;
2120  mOUT.p = NULL; mOUT.n = 0;
2121 
2122 
2123  INITIALIZE_GEOSPATIAL(node)
2124  COMPILE_GEOSYSTEM(node)
2125  MOVE_TO_ORIGIN(node)
2126 
2127 
2128  /* keep the output values of this process */
2129  FREE_IF_NZ(node->__movedValue.p);
2130  node->__movedValue.p = mOUT.p;
2131  node->__movedValue.n = mOUT.n;
2132 
2133  FREE_IF_NZ(gdCoords.p);
2134  MARK_NODE_COMPILED
2135 
2136  /* events */
2137  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoPositionInterpolator, metadata)) */
2138 }
2139 
2140 /* PositionInterpolator, ColorInterpolator, GeoPositionInterpolator */
2141 /* Called during the "events_processed" section of the event loop, */
2142 /* so this is called ONLY when there is something required to do, thus */
2143 /* there is no need to look at whether it is active or not */
2144 
2145 /* GeoPositionInterpolator == PositionIterpolator but with geovalue_changed and coordinate conversions */
2146 void do_GeoPositionInterpolator (void *innode) {
2147  struct X3D_GeoPositionInterpolator *node;
2148  int kin, kvin, counter, tmp;
2149  struct SFVec3d *kVs;
2150  /* struct SFColor *kVs */
2151 
2152  if (!innode) return;
2153  node = (struct X3D_GeoPositionInterpolator *) innode;
2154 
2155  if (NODE_NEEDS_COMPILING) compile_GeoPositionInterpolator(node);
2156  kvin = node->__movedValue.n;
2157  kVs = node->__movedValue.p;
2158  kin = node->key.n;
2159  MARK_EVENT (innode, offsetof (struct X3D_GeoPositionInterpolator, value_changed));
2160  MARK_EVENT (innode, offsetof (struct X3D_GeoPositionInterpolator, geovalue_changed));
2161 
2162  /* did the key or keyValue change? */
2163  if (node->__oldKeyValuePtr.p != node->keyValue.p) {
2164  MARK_EVENT (innode, offsetof (struct X3D_GeoPositionInterpolator, keyValue));
2165  node->__oldKeyValuePtr.p= node->keyValue.p;
2166  }
2167  if (node->__oldKeyPtr.p != node->key.p) {
2168  MARK_EVENT (innode, offsetof (struct X3D_GeoPositionInterpolator, key));
2169  node->__oldKeyPtr.p = node->key.p;
2170  }
2171 
2172 
2173  #ifdef SEVERBOSE
2174  printf("do_GeoPos: Position/Color interp, node %u kin %d kvin %d set_fraction %f\n",
2175  node, kin, kvin, node->set_fraction);
2176  #endif
2177 
2178  /* make sure we have the keys and keyValues */
2179  if ((kvin == 0) || (kin == 0)) {
2180  node->value_changed.c[0] = (float) 0.0;
2181  node->value_changed.c[1] = (float) 0.0;
2182  node->value_changed.c[2] = (float) 0.0;
2183  node->geovalue_changed.c[0] = 0.0;
2184  node->geovalue_changed.c[1] = 0.0;
2185  node->geovalue_changed.c[2] = 0.0;
2186  return;
2187  }
2188 
2189  if (kin>kvin) kin=kvin; /* means we don't use whole of keyValue, but... */
2190 
2191  /* set_fraction less than or greater than keys */
2192  if (node->set_fraction <= ((node->key).p[0])) {
2193  memcpy ((void *)&node->geovalue_changed, (void *)&kVs[0], sizeof (struct SFVec3d));
2194  } else if (node->set_fraction >= node->key.p[kin-1]) {
2195  memcpy ((void *)&node->geovalue_changed, (void *)&kVs[kvin-1], sizeof (struct SFVec3d));
2196  } else {
2197  /* have to go through and find the key before */
2198  counter = find_key(kin,((float)(node->set_fraction)),node->key.p);
2199  for (tmp=0; tmp<3; tmp++) {
2200  node->geovalue_changed.c[tmp] =
2201  (node->set_fraction - node->key.p[counter-1]) /
2202  (node->key.p[counter] - node->key.p[counter-1]) *
2203  (kVs[counter].c[tmp] - kVs[counter-1].c[tmp]) + kVs[counter-1].c[tmp];
2204  }
2205  }
2206 
2207  /* convert this back into the requested spatial format */
2208  CONVERT_BACK_TO_GD_OR_UTM(node->geovalue_changed)
2209 
2210  /* set the (float) value_changed, as well */
2211  for (tmp=0;tmp<3;tmp++) node->value_changed.c[tmp] = (float)node->geovalue_changed.c[tmp];
2212 
2213  #ifdef SEVERBOSE
2214  printf ("Pos/Col, new value (%f %f %f)\n",
2215  node->value_changed.c[0],node->value_changed.c[1],node->value_changed.c[2]);
2216  printf ("geovalue_changed %lf %lf %lf\n",node->geovalue_changed.c[0], node->geovalue_changed.c[1], node->geovalue_changed.c[2]);
2217  #endif
2218 }
2219 
2220 /************************************************************************/
2221 /* GeoProximitySensor */
2222 /************************************************************************/
2223 
2224 void compile_GeoProximitySensor (struct X3D_GeoProximitySensor * node) {
2225  MF_SF_TEMPS
2226 
2227  #ifdef VERBOSE
2228  printf ("compiling GeoProximitySensor\n");
2229  #endif
2230 
2231  /* work out the position */
2232  INITIALIZE_GEOSPATIAL(node)
2233  COMPILE_GEOSYSTEM(node)
2234  INIT_MF_FROM_SF(node, geoCenter)
2235  MOVE_TO_ORIGIN(node)
2236  COPY_MF_TO_SF(node, __movedCoords)
2237 
2238  /* work out the local orientation */
2239  GeoOrient(node->geoOrigin, &node->__geoSystem, &gdCoords.p[0], &node->__localOrient);
2240  #ifdef VERBOSE
2241  printf ("compile_GeoProximitySensor, orig coords %lf %lf %lf, moved %lf %lf %lf\n", node->geoCenter.c[0], node->geoCenter.c[1], node->geoCenter.c[2], node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
2242  printf (" rotation is %lf %lf %lf %lf\n",
2243  node->__localOrient.c[0],
2244  node->__localOrient.c[1],
2245  node->__localOrient.c[2],
2246  node->__localOrient.c[3]);
2247  #endif
2248 
2249  MARK_NODE_COMPILED
2250  FREE_MF_SF_TEMPS
2251 
2252  MARK_SFVEC3D_INOUT_EVENT(node->geoCenter, node->__oldGeoCenter,offsetof (struct X3D_GeoProximitySensor, geoCenter))
2253  MARK_SFVEC3F_INOUT_EVENT(node->size, node->__oldSize,offsetof (struct X3D_GeoProximitySensor, size))
2254 
2255  /* events */
2256  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoProximitySensor, metadata)) */
2257 
2258 
2259  #ifdef VERBOSE
2260  printf ("compiled GeoProximitySensor\n\n");
2261  #endif
2262 }
2263 
2264  //PROXIMITYSENSOR(GeoProximitySensor,__movedCoords,INITIALIZE_GEOSPATIAL(node),COMPILE_IF_REQUIRED)
2265 //#define PROXIMITYSENSOR(type,center,initializer1,initializer2)
2266 void proximity_GeoProximitySensor (struct X3D_GeoProximitySensor *node) {
2267  /* Viewer pos = t_r2 */
2268  double cx,cy,cz;
2269  double len;
2270  struct point_XYZ dr1r2;
2271  struct point_XYZ dr2r3;
2272  struct point_XYZ nor1,nor2;
2273  struct point_XYZ ins;
2274  static struct point_XYZ yvec = {0,0.05,0};
2275  static struct point_XYZ zvec = {0,0,-0.05};
2276  static struct point_XYZ zpvec = {0,0,0.05};
2277  static struct point_XYZ orig = {0,0,0};
2278  struct point_XYZ t_zvec, t_yvec, t_orig, t_center;
2279  GLDOUBLE modelMatrix[16];
2280  GLDOUBLE projMatrix[16];
2281  GLDOUBLE view2prox[16];
2282 
2283  if(!((node->enabled))) return;
2284  INITIALIZE_GEOSPATIAL(node)
2285  COMPILE_IF_REQUIRED
2286 
2287  /* printf (" vp %d geom %d light %d sens %d blend %d prox %d col %d\n",*/
2288  /* render_vp,render_geom,render_light,render_sensitive,render_blend,render_proximity,render_collision);*/
2289 
2290  /* transforms viewers coordinate space into sensors coordinate space.
2291  * this gives the orientation of the viewer relative to the sensor.
2292  */
2293  FW_GL_GETDOUBLEV(GL_MODELVIEW_MATRIX, modelMatrix);
2294  if(0){
2295  FW_GL_GETDOUBLEV(GL_PROJECTION_MATRIX, projMatrix);
2296  FW_GLU_UNPROJECT(orig.x,orig.y,orig.z,modelMatrix,projMatrix,viewport,
2297  &t_orig.x,&t_orig.y,&t_orig.z);
2298  FW_GLU_UNPROJECT(zvec.x,zvec.y,zvec.z,modelMatrix,projMatrix,viewport,
2299  &t_zvec.x,&t_zvec.y,&t_zvec.z);
2300  FW_GLU_UNPROJECT(yvec.x,yvec.y,yvec.z,modelMatrix,projMatrix,viewport,
2301  &t_yvec.x,&t_yvec.y,&t_yvec.z);
2302  VECDIFF(t_zvec, t_orig, dr1r2); /* Z axis */
2303  VECDIFF(t_yvec, t_orig, dr2r3); /* Y axis */
2304 
2305  }
2306  matinverseAFFINE(view2prox,modelMatrix);
2307  if(1){
2308  //feature-AFFINE_GLU_UNPROJECT
2309  transform(&t_orig,&orig,view2prox);
2310  transform(&zvec,&zvec,view2prox);
2311  transform(&yvec,&yvec,view2prox);
2312  VECDIFF(zvec, t_orig, dr1r2);
2313  VECDIFF(yvec, t_orig, dr2r3);
2314  }
2315  transform(&t_center,&orig, view2prox);
2316 
2317 
2318  /*printf ("\n");
2319  printf ("unprojected, t_orig (0,0,0) %lf %lf %lf\n",t_orig.x, t_orig.y, t_orig.z);
2320  printf ("unprojected, t_yvec (0,0.05,0) %lf %lf %lf\n",t_yvec.x, t_yvec.y, t_yvec.z);
2321  printf ("unprojected, t_zvec (0,0,-0.05) %lf %lf %lf\n",t_zvec.x, t_zvec.y, t_zvec.z);
2322  */
2323  cx = t_center.x - ((node->__movedCoords ).c[0]);
2324  cy = t_center.y - ((node->__movedCoords ).c[1]);
2325  cz = t_center.z - ((node->__movedCoords ).c[2]);
2326 
2327  if(((node->size).c[0]) == 0 || ((node->size).c[1]) == 0 || ((node->size).c[2]) == 0) return;
2328 
2329  if(fabs(cx) > ((node->size).c[0])/2 ||
2330  fabs(cy) > ((node->size).c[1])/2 ||
2331  fabs(cz) > ((node->size).c[2])/2) return;
2332  /* printf ("within (Geo)ProximitySensor\n"); */
2333 
2334  /* Ok, we now have to compute... */
2335  (node->__hit) /*cget*/ = 1;
2336 
2337  /* Position */
2338  ((node->__t1).c[0]) = (float)t_center.x;
2339  ((node->__t1).c[1]) = (float)t_center.y;
2340  ((node->__t1).c[2]) = (float)t_center.z;
2341 
2342 
2343  /* printf (" dr1r2 %lf %lf %lf\n",dr1r2.x, dr1r2.y, dr1r2.z);
2344  printf (" dr2r3 %lf %lf %lf\n",dr2r3.x, dr2r3.y, dr2r3.z);
2345  */
2346 
2347  len = sqrt(VECSQ(dr1r2)); VECSCALE(dr1r2,1/len);
2348  len = sqrt(VECSQ(dr2r3)); VECSCALE(dr2r3,1/len);
2349 
2350  /* printf ("scaled dr1r2 %lf %lf %lf\n",dr1r2.x, dr1r2.y, dr1r2.z);
2351  printf ("scaled dr2r3 %lf %lf %lf\n",dr2r3.x, dr2r3.y, dr2r3.z);
2352  */
2353 
2354  /*
2355  printf("PROX_INT: (%f %f %f) (%f %f %f) (%f %f %f)\n (%f %f %f) (%f %f %f)\n",
2356  t_orig.x, t_orig.y, t_orig.z,
2357  t_zvec.x, t_zvec.y, t_zvec.z,
2358  t_yvec.x, t_yvec.y, t_yvec.z,
2359  dr1r2.x, dr1r2.y, dr1r2.z,
2360  dr2r3.x, dr2r3.y, dr2r3.z
2361  );
2362  */
2363 
2364  if(fabs(VECPT(dr1r2, dr2r3)) > 0.001) {
2365  printf ("Sorry, can't handle unevenly scaled ProximitySensors yet :("
2366  "dp: %f v: (%f %f %f) (%f %f %f)\n", VECPT(dr1r2, dr2r3),
2367  dr1r2.x,dr1r2.y,dr1r2.z,
2368  dr2r3.x,dr2r3.y,dr2r3.z
2369  );
2370  return;
2371  }
2372 
2373 
2374  if(APPROX(dr1r2.z,1.0)) {
2375  /* rotation */
2376  ((node->__t2).c[0]) = (float) 0;
2377  ((node->__t2).c[1]) = (float) 0;
2378  ((node->__t2).c[2]) = (float) 1;
2379  ((node->__t2).c[3]) = (float) atan2(-dr2r3.x,dr2r3.y);
2380  } else if(APPROX(dr2r3.y,1.0)) {
2381  /* rotation */
2382  ((node->__t2).c[0]) = (float) 0;
2383  ((node->__t2).c[1]) = (float) 1;
2384  ((node->__t2).c[2]) = (float) 0;
2385  ((node->__t2).c[3]) = (float) atan2(dr1r2.x,dr1r2.z);
2386  } else {
2387  /* Get the normal vectors of the possible rotation planes */
2388  nor1 = dr1r2;
2389  nor1.z -= 1.0;
2390  nor2 = dr2r3;
2391  nor2.y -= 1.0;
2392 
2393  /* Now, the intersection of the planes, obviously cp */
2394  VECCP(nor1,nor2,ins);
2395 
2396  len = sqrt(VECSQ(ins)); VECSCALE(ins,1/len);
2397 
2398  /* the angle */
2399  VECCP(dr1r2,ins, nor1);
2400  VECCP(zpvec, ins, nor2);
2401  len = sqrt(VECSQ(nor1)); VECSCALE(nor1,1/len);
2402  len = sqrt(VECSQ(nor2)); VECSCALE(nor2,1/len);
2403  VECCP(nor1,nor2,ins);
2404 
2405  ((node->__t2).c[3]) = (float) -atan2(sqrt(VECSQ(ins)), VECPT(nor1,nor2));
2406 
2407  /* rotation - should normalize sometime... */
2408  ((node->__t2).c[0]) = (float) ins.x;
2409  ((node->__t2).c[1]) = (float) ins.y;
2410  ((node->__t2).c[2]) = (float) ins.z;
2411  }
2412  /*
2413  printf("NORS: (%f %f %f) (%f %f %f) (%f %f %f)\n",
2414  nor1.x, nor1.y, nor1.z,
2415  nor2.x, nor2.y, nor2.z,
2416  ins.x, ins.y, ins.z
2417  );
2418  */
2419 }
2420 
2421 
2422 /* GeoProximitySensor code for ClockTick */
2423 void do_GeoProximitySensorTick( void *ptr) {
2424  struct X3D_GeoProximitySensor *node = (struct X3D_GeoProximitySensor *)ptr;
2425 
2426  /* if not enabled, do nothing */
2427  if (!node) return;
2428  if (node->__oldEnabled != node->enabled) {
2429  node->__oldEnabled = node->enabled;
2430  MARK_EVENT(X3D_NODE(node),offsetof (struct X3D_GeoProximitySensor, enabled));
2431  }
2432  if (!node->enabled) return;
2433 
2434  /* isOver state */
2435  /* did we get a signal? */
2436  if (node->__hit) {
2437  if (!node->isActive) {
2438  #ifdef SEVERBOSE
2439  printf ("PROX - initial defaults\n");
2440  #endif
2441 
2442  node->isActive = 1;
2443  node->enterTime = TickTime();
2444  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, isActive));
2445  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, enterTime));
2446 
2447  }
2448 
2449  /* now, has anything changed? */
2450  if (memcmp ((void *) &node->position_changed,(void *) &node->__t1,sizeof(struct SFColor))) {
2451  #ifdef SEVERBOSE
2452  printf ("PROX - position changed!!! \n");
2453  #endif
2454 
2455  memcpy ((void *) &node->position_changed,
2456  (void *) &node->__t1,sizeof(struct SFColor));
2457  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, position_changed));
2458 
2459  #ifdef VERBOSE
2460  printf ("do_GeoProximitySensorTick, position changed; it now is %lf %lf %lf\n",node->position_changed.c[0],
2461  node->position_changed.c[1], node->position_changed.c[2]);
2462  printf ("nearPlane is %lf\n",Viewer.nearPlane);
2463 
2464  #endif
2465 
2466  /* possibly we have to convert this from GCC to GDC, and maybe even then to UTM */
2467 
2468  /* prep the geoCoord changed; first, get the position. Right now, we use the
2469  Viewer position, as it is more accurate (not clipped by the nearPlane) than
2470  the position_changed field */
2471 
2472  node->geoCoord_changed.c[0] = (double) node->position_changed.c[0];
2473  node->geoCoord_changed.c[1] = (double) node->position_changed.c[1];
2474  node->geoCoord_changed.c[2] = (double) node->position_changed.c[2];
2475 
2476  /* then add in the nearPlane, as the way we get the position is via a clipped frustum */
2477  /* if we get this via the position_changed field, we have to:
2478  node->geoCoord_changed.c[2] += Viewer.nearPlane;
2479  */
2480  node->geoCoord_changed.c[2] += Viewer()->nearPlane;
2481  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, geoCoord_changed));
2482 
2483  #ifdef VERBOSE
2484  printf ("\ngeoCoord_changed as a GCC, %lf %lf %lf\n",
2485  node->geoCoord_changed.c[0],
2486  node->geoCoord_changed.c[1],
2487  node->geoCoord_changed.c[2]);
2488  #endif
2489 
2490  CONVERT_BACK_TO_GD_OR_UTM(node->geoCoord_changed)
2491  }
2492  if (memcmp ((void *) &node->orientation_changed, (void *) &node->__t2,sizeof(struct SFRotation))) {
2493  #ifdef SEVERBOSE
2494  printf ("PROX - orientation changed!!!\n ");
2495  #endif
2496 
2497  memcpy ((void *) &node->orientation_changed,
2498  (void *) &node->__t2,sizeof(struct SFRotation));
2499  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, orientation_changed));
2500  }
2501  } else {
2502  if (node->isActive) {
2503  #ifdef SEVERBOSE
2504  printf ("PROX - stopping\n");
2505  #endif
2506 
2507  node->isActive = 0;
2508  node->exitTime = TickTime();
2509  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, isActive));
2510 
2511  MARK_EVENT (ptr, offsetof(struct X3D_GeoProximitySensor, exitTime));
2512  }
2513  }
2514  node->__hit=FALSE;
2515 }
2516 
2517 
2518 /************************************************************************/
2519 /* GeoTouchSensor */
2520 /************************************************************************/
2521 
2522 void compile_GeoTouchSensor (struct X3D_GeoTouchSensor * node) {
2523  #ifdef VERBOSE
2524  printf ("compiling GeoTouchSensor\n");
2525  #endif
2526 
2527  INITIALIZE_GEOSPATIAL(node)
2528  COMPILE_GEOSYSTEM(node)
2529  MARK_NODE_COMPILED
2530 
2531  /* events */
2532  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoTouchSensor, metadata)) */
2533 
2534 }
2535 
2536 void do_GeoTouchSensor ( void *ptr, int ev, int but1, int over) {
2537 
2538 
2539  struct X3D_GeoTouchSensor *node = (struct X3D_GeoTouchSensor *)ptr;
2540  struct point_XYZ normalval; /* different structures for normalization calls */
2541  ttglobal tg;
2542  COMPILE_IF_REQUIRED
2543 
2544  #ifdef SENSVERBOSE
2545  printf ("%lf: TS ",TickTime());
2546  if (ev==ButtonPress) printf ("ButtonPress ");
2547  else if (ev==ButtonRelease) printf ("ButtonRelease ");
2548  else if (ev==KeyPress) printf ("KeyPress ");
2549  else if (ev==KeyRelease) printf ("KeyRelease ");
2550  else if (ev==MotionNotify) printf ("%lf MotionNotify ");
2551  else printf ("ev %d ",ev);
2552 
2553  if (but1) printf ("but1 TRUE "); else printf ("but1 FALSE ");
2554  if (over) printf ("over TRUE "); else printf ("over FALSE ");
2555  printf ("\n");
2556  #endif
2557 
2558  /* if not enabled, do nothing */
2559  if (!node) return;
2560  if (node->__oldEnabled != node->enabled) {
2561  node->__oldEnabled = node->enabled;
2562  MARK_EVENT(X3D_NODE(node),offsetof (struct X3D_GeoTouchSensor, enabled));
2563  }
2564  if (!node->enabled) return;
2565  tg = gglobal();
2566  /* isOver state */
2567  if ((ev == overMark) && (over != node->isOver)) {
2568  #ifdef SENSVERBOSE
2569  printf ("TS %u, isOver changed %d\n",node, over);
2570  #endif
2571  node->isOver = over;
2572  MARK_EVENT (ptr, offsetof (struct X3D_GeoTouchSensor, isOver));
2573  }
2574 
2575  /* active */
2576  /* button presses */
2577  if (ev == ButtonPress) {
2578  node->isActive=1;
2579  MARK_EVENT (ptr, offsetof (struct X3D_GeoTouchSensor, isActive));
2580  #ifdef SENSVERBOSE
2581  printf ("touchSens %u, butPress\n",node);
2582  #endif
2583 
2584  node->touchTime = TickTime();
2585  MARK_EVENT(ptr, offsetof (struct X3D_GeoTouchSensor, touchTime));
2586 
2587  } else if (ev == ButtonRelease) {
2588  #ifdef SENSVERBOSE
2589  printf ("touchSens %u, butRelease\n",node);
2590  #endif
2591  node->isActive=0;
2592  MARK_EVENT (ptr, offsetof (struct X3D_GeoTouchSensor, isActive));
2593  }
2594 
2595  /* hitPoint and hitNormal */
2596  /* save the current hitPoint for determining if this changes between runs */
2597  memcpy ((void *) &node->_oldhitPoint, (void *) &tg->RenderFuncs.ray_save_posn,sizeof(struct SFColor));
2598 
2599  /* did the hitPoint change between runs? */
2600  if ((APPROX(node->_oldhitPoint.c[0],node->hitPoint_changed.c[0])!= TRUE) ||
2601  (APPROX(node->_oldhitPoint.c[1],node->hitPoint_changed.c[1])!= TRUE) ||
2602  (APPROX(node->_oldhitPoint.c[2],node->hitPoint_changed.c[2])!= TRUE)) {
2603 
2604  #ifdef SENSVERBOSE
2605  printf ("GeoTouchSens, hitPoint changed: %f %f %f\n",node->hitPoint_changed.c[0],
2606  node->hitPoint_changed.c[1], node->hitPoint_changed.c[2]);
2607  #endif
2608 
2609  memcpy ((void *) &node->hitPoint_changed, (void *) &node->_oldhitPoint, sizeof(struct SFColor));
2610  MARK_EVENT(ptr, offsetof (struct X3D_GeoTouchSensor, hitPoint_changed));
2611 
2612  /* convert this back into the requested GeoSpatial format... */
2613  node->hitGeoCoord_changed.c[0] = (double) node->hitPoint_changed.c[0];
2614  node->hitGeoCoord_changed.c[1] = (double) node->hitPoint_changed.c[1];
2615  node->hitGeoCoord_changed.c[2] = (double) node->hitPoint_changed.c[2];
2616 
2617  /* then add in the nearPlane, as the way we get the position is via a clipped frustum */
2618  /* if we get this via the position_changed field, we have to:
2619  node->hitGeoCoord_changed.c[2] += nearPlane;
2620  */
2621  node->hitGeoCoord_changed.c[2] += Viewer()->nearPlane;
2622  MARK_EVENT (ptr, offsetof(struct X3D_GeoTouchSensor, hitGeoCoord_changed));
2623 
2624  #ifdef SENSVERBOSE
2625  printf ("\nhitGeoCoord_changed as a GCC, %lf %lf %lf\n",
2626  node->hitGeoCoord_changed.c[0],
2627  node->hitGeoCoord_changed.c[1],
2628  node->hitGeoCoord_changed.c[2]);
2629  #endif
2630 
2631  CONVERT_BACK_TO_GD_OR_UTM(node->hitGeoCoord_changed)
2632  }
2633 
2634  /* have to normalize normal; change it from SFColor to struct point_XYZ. */
2635  normalval.x = tg->RenderFuncs.hyp_save_norm[0];
2636  normalval.y = tg->RenderFuncs.hyp_save_norm[1];
2637  normalval.z = tg->RenderFuncs.hyp_save_norm[2];
2638  normalize_vector(&normalval);
2639  node->_oldhitNormal.c[0] = (float) normalval.x;
2640  node->_oldhitNormal.c[1] = (float) normalval.y;
2641  node->_oldhitNormal.c[2] = (float) normalval.z;
2642 
2643  /* did the hitNormal change between runs? */
2644  MARK_SFVEC3F_INOUT_EVENT(node->hitNormal_changed,node->_oldhitNormal,offsetof (struct X3D_GeoTouchSensor, hitNormal_changed))
2645 }
2646 
2647 
2648 
2649 /************************************************************************/
2650 /* GeoViewpoint */
2651 /************************************************************************/
2652 
2653 void compile_GeoViewpoint (struct X3D_GeoViewpoint * node) {
2654  struct SFVec4d localOrient;
2655  struct SFVec4d orient;
2656  int i;
2657  Quaternion localQuat;
2658  Quaternion relQuat;
2659  Quaternion combQuat;
2660  MF_SF_TEMPS
2661 
2662  #ifdef VERBOSE
2663  printf ("compileViewpoint is %u, its geoOrigin is %u \n",node, node->geoOrigin);
2664  if (node->geoOrigin!=NULL) printf ("type %s\n",stringNodeType(X3D_GEOORIGIN(node->geoOrigin)->_nodeType));
2665  #endif
2666 
2667 
2668  /* did any of the "set_" inputOnly fields get set? if not, just use the non-set fields */
2669  USE_SET_SFVEC3D_IF_CHANGED(set_position,position)
2670  USE_SET_SFROTATION_IF_CHANGED(set_orientation,orientation)
2671 
2672  /* work out the position */
2673  INITIALIZE_GEOSPATIAL(node)
2674  COMPILE_GEOSYSTEM(node)
2675  INIT_MF_FROM_SF(node, position)
2676  MOVE_TO_ORIGIN(node)
2677  COPY_MF_TO_SF(node, __movedPosition)
2678 
2679  /* work out the local orientation and copy doubles to floats */
2680  GeoOrient(node->geoOrigin, &node->__geoSystem, &gdCoords.p[0], &localOrient);
2681 
2682  /* Quaternize the local Geospatial quaternion, and the specified rotation from the GeoViewpoint orientation field */
2683  vrmlrot_to_quaternion (&localQuat, localOrient.c[0], localOrient.c[1], localOrient.c[2], localOrient.c[3]);
2684  vrmlrot_to_quaternion (&relQuat, node->orientation.c[0], node->orientation.c[1], node->orientation.c[2], node->orientation.c[3]);
2685 
2686  /* add these together */
2687  quaternion_add (&combQuat, &relQuat, &localQuat);
2688 
2689  /* get the rotation; 2 steps to convert doubles to floats;
2690  should be quaternion_to_vrmlrot(&combQuat, &node->__movedOrientation.c[0]... */
2691  quaternion_to_vrmlrot(&combQuat, &orient.c[0], &orient.c[1], &orient.c[2], &orient.c[3]);
2692  for (i=0; i<4; i++) node->__movedOrientation.c[i] = (float) orient.c[i];
2693 
2694  #ifdef VERBOSE
2695  printf ("compile_GeoViewpoint, final position %lf %lf %lf\n",node->__movedPosition.c[0],
2696  node->__movedPosition.c[1], node->__movedPosition.c[2]);
2697 
2698  printf ("compile_GeoViewpoint, getLocalOrientation %lf %lf %lf %lf\n",localOrient.c[0],
2699  localOrient.c[1], localOrient.c[2], localOrient.c[3]);
2700  printf ("compile_GeoViewpoint, initial orientation: %lf %lf %lf %lf\n",node->orientation.c[0],
2701  node->orientation.c[1], node->orientation.c[2], node->orientation.c[3]);
2702  printf ("compile_GeoViewpoint, final rotation %lf %lf %lf %lf\n",node->__movedOrientation.c[0],
2703  node->__movedOrientation.c[1], node->__movedOrientation.c[2], node->__movedOrientation.c[3]);
2704  printf ("compile_GeoViewpoint, elevation from the WGS84 ellipsoid is %lf\n",gdCoords.p[0].c[2]);
2705  #endif
2706 
2707  MARK_NODE_COMPILED
2708  FREE_MF_SF_TEMPS
2709 
2710  /* events */
2711  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoViewpoint, metadata)) */
2712  MARK_SFFLOAT_INOUT_EVENT(node->fieldOfView, node->__oldFieldOfView, offsetof (struct X3D_GeoViewpoint, fieldOfView))
2713  MARK_SFBOOL_INOUT_EVENT(node->headlight, node->__oldHeadlight, offsetof (struct X3D_GeoViewpoint, headlight))
2714  MARK_SFBOOL_INOUT_EVENT(node->jump, node->__oldJump, offsetof (struct X3D_GeoViewpoint, jump))
2715  /*
2716  //dug9 may 2015 I'm not sure what the __old stuff was for (H: debugging) but
2717  //shallow copying a string pointer -or MFString or SFString- makes it hard to generically free during exit
2718  //see opengl_utils.c in cbFreeMallocedBuiltinField
2719  MARK_SFSTRING_INOUT_EVENT(node->description,node->__oldSFString, offsetof(struct X3D_GeoViewpoint, description))
2720  MARK_MFSTRING_INOUT_EVENT(node->navType,node->__oldMFString, offsetof(struct X3D_GeoViewpoint, navType))
2721  */
2722  #ifdef VERBOSE
2723  printf ("compiled GeoViewpoint\n\n");
2724  #endif
2725 }
2726 struct X3D_Node *getActiveLayerBoundViewpoint();
2727 void prep_GeoViewpoint (struct X3D_GeoViewpoint *node) {
2728  double a1;
2729  GLint viewPort[10];
2730  if (!renderstate()->render_vp) return;
2731 
2732  if((struct X3D_Node*)node == getActiveLayerBoundViewpoint() && !node->_donethispass){
2733  node->_donethispass = 1; //if the vp id DEF/USED multiple places in the scengraph,
2734 
2735  INITIALIZE_GEOSPATIAL(node)
2736 
2737  /* printf ("RVP, node %d ib %d sb %d gepvp\n",node,node->isBound,node->set_bind);
2738  printf ("VP stack %d tos %d\n",viewpoint_tos, viewpoint_stack[viewpoint_tos]);
2739  */
2740 
2741  /* check the set_bind eventin to see if it is TRUE or FALSE */
2742  /* code to perform binding is now in set_viewpoint. */
2743 
2744  COMPILE_IF_REQUIRED
2745 
2746  #ifdef VERBOSE
2747  printf ("prep_GeoViewpoint called\n");
2748  #endif
2749 
2750  /* perform GeoViewpoint translations */
2751  FW_GL_ROTATE_RADIANS(-node->__movedOrientation.c[3],node->__movedOrientation.c[0],node->__movedOrientation.c[1],
2752  node->__movedOrientation.c[2]);
2753 
2754  FW_GL_TRANSLATE_D(-node->__movedPosition.c[0],-node->__movedPosition.c[1],-node->__movedPosition.c[2]);
2755 
2756  /* we have a new currentPosInModel now... */
2757  /* printf ("currentPosInModel was %lf %lf %lf\n", Viewer.currentPosInModel.x, Viewer.currentPosInModel.y, Viewer.currentPosInModel.z); */
2758 
2759  /* the AntiPos has been applied in the trans and rots above, so we do not need to do it here */
2760  getCurrentPosInModel(FALSE);
2761 
2762 
2763  /* now, lets work on the GeoViewpoint fieldOfView */
2764  FW_GL_GETINTEGERV(GL_VIEWPORT, viewPort);
2765  if(viewPort[2] > viewPort[3]) {
2766  a1=0;
2767  Viewer()->fieldofview = node->fieldOfView/3.1415926536*180;
2768  } else {
2769  a1 = node->fieldOfView;
2770  a1 = atan2(sin(a1),viewPort[2]/((float)viewPort[3]) * cos(a1));
2771  Viewer()->fieldofview = a1/3.1415926536*180;
2772  }
2773 
2774  calculateViewingSpeed();
2775  #ifdef VERBOSE
2776  printf ("prep_GeoViewpoint, fieldOfView %f \n",node->fieldOfView);
2777  #endif
2778  }
2779 }
2780 
2781 /* GeoViewpoint speeds and avatar sizes are depenent on elevation above WGS_84. These are calculated here */
2782 /* this is called from the Viewer functions */
2783 void calculateViewingSpeed() {
2784  struct SFVec3d gcCoords;
2785  struct SFVec3d gdCoords;
2786 
2787  /* the current position is the GC coordinate */
2788  gcCoords.c[0]= Viewer()->currentPosInModel.x;
2789  gcCoords.c[1] = Viewer()->currentPosInModel.y;
2790  gcCoords.c[2] = Viewer()->currentPosInModel.z;
2791 
2792  #ifdef VERBOSE
2793  printf ("calculateViewingSpeed, currentPosInModel %lf %lf %lf\n", gcCoords.c[0], gcCoords.c[1], gcCoords.c[2]);
2794  #endif
2795 
2796  if (Viewer()->GeoSpatialNode != NULL) {
2797  /* do we have a valid __geoSystem?? */
2798  INITIALIZE_GEOSPATIAL(Viewer()->GeoSpatialNode)
2799 
2800 
2801 /*
2802  COMPILE_IF_REQUIRED
2803 
2804 */
2805 
2806 
2807 
2808 
2809  if (Viewer()->GeoSpatialNode->__geoSystem.n>0) {
2810  /* is the __geoSystem NOT gc coords? */
2811  /* printf ("have a GeoSpatial viewpoint, currently %d\n",Viewer.GeoSpatialNode->__geoSystem.p[0]); */
2812  if (Viewer()->GeoSpatialNode->__geoSystem.p[0] != GEOSP_GC) {
2813 
2814 /*
2815  retractOrigin((struct X3D_GeoOrigin *)Viewer.GeoSpatialNode->geoOrigin, &gcCoords);
2816 */
2817 
2818  #ifdef VERBOSE
2819  printf ("\n");
2820  printf ("for GeoViewpoint :%s:\n",Viewer.GeoSpatialNode->description->strptr);
2821  printf ("calculateViewingSpeed, currentPosInModel: %lf %lf %lf\n", gcCoords.c[0], gcCoords.c[1], gcCoords.c[2]);
2822  #endif
2823 
2824  /* convert from local (gc) to gd coordinates, using WGS84 ellipsoid */
2825  gccToGdc (&gcCoords, &gdCoords);
2826 
2827  #ifdef VERBOSE
2828  printf ("speed is calculated from geodetic height %lf %lf %lf\n",gdCoords.c[0], gdCoords.c[1], gdCoords.c[2]);
2829  #endif
2830 
2831  /* speed is dependent on elevation above WGS84 ellipsoid */
2832  Viewer()->speed = fabs(sqrt(gcCoords.c[0]*gcCoords.c[0] + gcCoords.c[1]*gcCoords.c[1] + gcCoords.c[2]*gcCoords.c[2])
2833  -GEOSP_WE_A) * Viewer()->GeoSpatialNode->speedFactor;
2834  if (Viewer()->speed < 1.0) Viewer()->speed=1.0;
2835 
2836  #ifdef VERBOSE
2837  printf ("height above center %f WGS84 ellipsoid is %lf\n",Viewer.speed,GEOSP_WE_A);
2838  #endif
2839 /*
2840  Viewer.speed = fabs(Viewer.speed * Viewer.GeoSpatialNode->speedFactor);
2841  if (Viewer.speed < Viewer.GeoSpatialNode->speedFactor) Viewer.speed = Viewer.GeoSpatialNode->speedFactor;
2842 */
2843 
2844  /* set the navigation info - use the GeoVRML algorithms */
2845  set_naviWidthHeightStep(
2846  Viewer()->speed*0.25,
2847  Viewer()->speed*1.6,
2848  Viewer()->speed*0.25);
2849  }
2850  }
2851  }
2852 }
2853 
2854 static void calculateExamineModeDistance(void) {
2855 /*
2856  printf ("bind_GeoViewpoint - calculateExamineModeDistance\n");
2857 */
2858 Viewer()->doExamineModeDistanceCalculations = TRUE;
2859 
2860 }
2861 
2862 void bind_GeoViewpoint (struct X3D_GeoViewpoint *node) {
2863  X3D_Viewer *viewer;
2864  Quaternion q_i;
2865 
2866  /* did bind_node tell us we could bind this guy? */
2867  if (!(node->isBound)) return;
2868 
2869  viewer = ViewerByLayerId(node->_layerId);
2870 
2871  INITIALIZE_GEOSPATIAL(node)
2872  COMPILE_IF_REQUIRED
2873 
2874  /* set Viewer position and orientation */
2875 
2876  #ifdef VERBOSE
2877  printf ("bind_GeoViewpoint, setting Viewer to %lf %lf %lf orient %f %f %f %f\n",node->__movedPosition.c[0],node->__movedPosition.c[1],
2878  node->__movedPosition.c[2],node->orientation.c[0],node->orientation.c[1],node->orientation.c[2],
2879  node->orientation.c[3]);
2880  printf (" node %u fieldOfView %f\n",node,node->fieldOfView);
2881  #endif
2882 
2883  viewer->GeoSpatialNode = node;
2884 
2885  viewer->Pos.x = node->__movedPosition.c[0];
2886  viewer->Pos.y = node->__movedPosition.c[1];
2887  viewer->Pos.z = node->__movedPosition.c[2];
2888  viewer->AntiPos.x = node->__movedPosition.c[0];
2889  viewer->AntiPos.y = node->__movedPosition.c[1];
2890  viewer->AntiPos.z = node->__movedPosition.c[2];
2891 
2892  /* printf ("bind_GeoViewpoint, pos %f %f %f antipos %f %f %f\n",Viewer.Pos.x, Viewer.Pos.y, Viewer.Pos.z, Viewer.AntiPos.x, Viewer.AntiPos.y, Viewer.AntiPos.z); */
2893 
2894  vrmlrot_to_quaternion (&viewer->Quat,node->__movedOrientation.c[0],
2895  node->__movedOrientation.c[1],node->__movedOrientation.c[2],node->__movedOrientation.c[3]);
2896 
2897  vrmlrot_to_quaternion (&q_i,node->__movedOrientation.c[0],
2898  node->__movedOrientation.c[1],node->__movedOrientation.c[2],node->__movedOrientation.c[3]);
2899  quaternion_inverse(&(viewer->AntiQuat),&q_i);
2900 
2901  resolve_pos();
2902 
2903  calculateViewingSpeed();
2904 
2905  calculateExamineModeDistance();
2906  setMenuStatusVP (node->description->strptr);
2907 
2908 }
2909 
2910 
2911 /************************************************************************/
2912 /* GeoTransform */
2913 /************************************************************************/
2914 
2915 void compile_GeoTransform (struct X3D_GeoTransform * node) {
2916  MF_SF_TEMPS
2917 
2918  #ifdef VERBOSE
2919  printf ("compiling GeoLocation\n");
2920  #endif
2921 
2922  /* work out the position */
2923  INITIALIZE_GEOSPATIAL(node)
2924  COMPILE_GEOSYSTEM(node)
2925  INIT_MF_FROM_SF(node, geoCenter)
2926  MOVE_TO_ORIGIN(node)
2927  COPY_MF_TO_SF(node, __movedCoords)
2928 
2929  /* work out the local orientation */
2930  GeoOrient(node->geoOrigin, &node->__geoSystem, &gdCoords.p[0], &node->__localOrient);
2931 
2932  MARK_SFVEC3D_INOUT_EVENT(node->geoCenter, node->__oldGeoCenter,offsetof (struct X3D_GeoTransform, geoCenter))
2933  MARK_MFNODE_INOUT_EVENT(node->children, node->__oldChildren, offsetof (struct X3D_GeoTransform, children))
2934 
2935 
2936  /* re-figure out which modifiers are actually in use */
2937  /* printf ("re-rendering for %d\n",node);*/
2938  node->__do_trans = verify_translate ((GLfloat *)node->translation.c);
2939  if (node->__do_trans) MARK_EVENT(X3D_NODE(node), offsetof (struct X3D_GeoTransform, translation));
2940 
2941  node->__do_scale = verify_scale ((GLfloat *)node->scale.c);
2942  if (node->__do_scale) MARK_EVENT(X3D_NODE(node), offsetof (struct X3D_GeoTransform, scale));
2943 
2944  node->__do_rotation = verify_rotate ((GLfloat *)node->rotation.c);
2945  if (node->__do_rotation) MARK_EVENT(X3D_NODE(node), offsetof (struct X3D_GeoTransform, rotation));
2946 
2947  node->__do_scaleO = verify_rotate ((GLfloat *)node->scaleOrientation.c);
2948  if (node->__do_scaleO) MARK_EVENT(X3D_NODE(node), offsetof (struct X3D_GeoTransform, scaleOrientation));
2949 
2950 
2951 
2952  #ifdef VERBOSE
2953  printf ("compile_GeoTransform, orig coords %lf %lf %lf, moved %lf %lf %lf\n", node->geoCoords.c[0], node->geoCoords.c[1], node->geoCoords.c[2], node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
2954  printf (" rotation is %lf %lf %lf %lf\n",
2955  node->__localOrient.c[0],
2956  node->__localOrient.c[1],
2957  node->__localOrient.c[2],
2958  node->__localOrient.c[3]);
2959  #endif
2960 
2961  REINITIALIZE_SORTED_NODES_FIELD(node->children,node->_sortedChildren);
2962  MARK_NODE_COMPILED
2963  FREE_MF_SF_TEMPS
2964 
2965  /* events */
2966  /* MARK_SFNODE_INOUT_EVENT(node->metadata, node->__oldmetadata, offsetof (struct X3D_GeoTransform, metadata)) */
2967 
2968 
2969  #ifdef VERBOSE
2970  printf ("compiled GeoTransform\n\n");
2971  #endif
2972 }
2973 
2974 
2975 /* do transforms, calculate the distance */
2976 void prep_GeoTransform (struct X3D_GeoTransform *node) {
2977 
2978  INITIALIZE_GEOSPATIAL(node)
2979  COMPILE_IF_REQUIRED
2980 
2981  /* rendering the viewpoint means doing the inverse transformations in reverse order (while poping stack),
2982  * so we do nothing here in that case -ncoder */
2983 
2984  /* printf ("prep_Transform, render_hier vp %d geom %d light %d sens %d blend %d prox %d col %d\n",
2985  render_vp,render_geom,render_light,render_sensitive,render_blend,render_proximity,render_collision); */
2986 
2987  /* do we have any geometry visible, and are we doing anything with geometry? */
2988  OCCLUSIONTEST
2989 
2990  if(!renderstate()->render_vp) {
2991  FW_GL_PUSH_MATRIX();
2992 
2993 
2994  /* TRANSLATION */
2995  if (node->__do_trans)
2996  FW_GL_TRANSLATE_F(node->translation.c[0],node->translation.c[1],node->translation.c[2]);
2997 
2998  /* GeoTransform TRANSLATION */
2999  FW_GL_TRANSLATE_D(node->__movedCoords.c[0], node->__movedCoords.c[1], node->__movedCoords.c[2]);
3000 
3001  //printf ("prep_GeoLoc trans to %lf %lf %lf\n",node->__movedCoords.c[0],node->__movedCoords.c[1],node->__movedCoords.c[2]);
3002 
3003 
3004  FW_GL_ROTATE_RADIANS(node->__localOrient.c[3], node->__localOrient.c[0],node->__localOrient.c[1],node->__localOrient.c[2]);
3005 
3006  /* ROTATION */
3007  if (node->__do_rotation) {
3008  FW_GL_ROTATE_RADIANS(node->rotation.c[3], node->rotation.c[0],node->rotation.c[1],node->rotation.c[2]);
3009  }
3010 
3011  /* SCALEORIENTATION */
3012  if (node->__do_scaleO) {
3013  FW_GL_ROTATE_RADIANS(node->scaleOrientation.c[3], node->scaleOrientation.c[0],
3014  node->scaleOrientation.c[1],node->scaleOrientation.c[2]);
3015  }
3016 
3017  /* SCALE */
3018  if (node->__do_scale)
3019  FW_GL_SCALE_F(node->scale.c[0],node->scale.c[1],node->scale.c[2]);
3020 
3021  /* REVERSE SCALE ORIENTATION */
3022  if (node->__do_scaleO)
3023  FW_GL_ROTATE_RADIANS(-node->scaleOrientation.c[3], node->scaleOrientation.c[0],
3024  node->scaleOrientation.c[1],node->scaleOrientation.c[2]);
3025 
3026  /* REVERSE CENTER */
3027  FW_GL_TRANSLATE_D(-node->__movedCoords.c[0], -node->__movedCoords.c[1], -node->__movedCoords.c[2]);
3028 
3029  RECORD_DISTANCE
3030  }
3031 }
3032 
3033 
3034 void fin_GeoTransform (struct X3D_GeoTransform *node) {
3035  INITIALIZE_GEOSPATIAL(node)
3036  COMPILE_IF_REQUIRED
3037  OCCLUSIONTEST
3038 
3039  if(!renderstate()->render_vp) {
3040  FW_GL_POP_MATRIX();
3041  } else {
3042  /*Rendering the viewpoint only means finding it, and calculating the reverse WorldView matrix.*/
3043  if((node->_renderFlags & VF_Viewpoint) == VF_Viewpoint) {
3044  FW_GL_TRANSLATE_D(((node->__movedCoords).c[0]),((node->__movedCoords).c[1]),((node->__movedCoords).c[2])
3045  );
3046  FW_GL_ROTATE_RADIANS(node->scaleOrientation.c[3],node->scaleOrientation.c[0],node->scaleOrientation.c[1],node->scaleOrientation.c[2]);
3047  FW_GL_SCALE_F((float)1.0/(((node->scale).c[0])),(float)1.0/(((node->scale).c[1])),(float)1.0/(((node->scale).c[2]))
3048  );
3049  FW_GL_ROTATE_RADIANS(-node->scaleOrientation.c[3],node->scaleOrientation.c[0],node->scaleOrientation.c[1],node->scaleOrientation.c[2]);
3050  FW_GL_ROTATE_RADIANS(-(((node->rotation).c[3])),((node->rotation).c[0]),((node->rotation).c[1]),((node->rotation).c[2])
3051  );
3052  FW_GL_TRANSLATE_D(-(((node->__movedCoords).c[0])),-(((node->__movedCoords).c[1])),-(((node->__movedCoords).c[2]))
3053  );
3054  FW_GL_TRANSLATE_F(-(((node->translation).c[0])),-(((node->translation).c[1])),-(((node->translation).c[2]))
3055  );
3056  }
3057  }
3058 }
3059 
3060 void child_GeoTransform (struct X3D_GeoTransform *node) {
3061  CHILDREN_COUNT
3062  //LOCAL_LIGHT_SAVE
3063  INITIALIZE_GEOSPATIAL(node)
3064  COMPILE_IF_REQUIRED
3065  OCCLUSIONTEST
3066  RETURN_FROM_CHILD_IF_NOT_FOR_ME
3067 
3068  /* any children at all? */
3069  if (nc==0) return;
3070 
3071  /* {
3072  int x;
3073  struct X3D_Node *xx;
3074 
3075  printf ("child_Transform, this %d \n",node);
3076  for (x=0; x<nc; x++) {
3077  xx = X3D_NODE(node->children.p[x]);
3078  printf (" ch %d type %s dist %f\n",node->children.p[x],stringNodeType(xx->_nodeType),xx->_dist);
3079  }
3080  } */
3081 
3082 
3083  /* do we have a local light for a child? */
3084  //LOCAL_LIGHT_CHILDREN(node->children);
3085  prep_sibAffectors((struct X3D_Node*)node,&node->__sibAffectors);
3086 
3087  /* now, just render the non-directionalLight children */
3088 
3089  /* printf ("Transform %d, flags %d, render_sensitive %d\n",
3090  node,node->_renderFlags,render_sensitive); */
3091 
3092  #ifdef CHILDVERBOSE
3093  printf ("transform - doing normalChildren\n");
3094  #endif
3095 
3096  normalChildren(node->children);
3097 
3098  #ifdef CHILDVERBOSE
3099  printf ("transform - done normalChildren\n");
3100  #endif
3101 
3102  //LOCAL_LIGHT_OFF
3103  fin_sibAffectors((struct X3D_Node*)node,&node->__sibAffectors);
3104 
3105 }
X3D_GeoPositionInterpolator
Definition: Structs.h:4885
X3D_GeoOrigin
Definition: Structs.h:4858
X3D_Node
Definition: Structs.h:2502
X3D_GeoLOD
Definition: Structs.h:4758
viewer
Definition: Viewer.h:196
X3D_GeoCoordinate
Definition: Structs.h:4695
X3D_GeoProximitySensor
Definition: Structs.h:4915
X3D_Proto
Definition: Structs.h:8567
X3D_PolyRep
Definition: Structs.h:62
X3D_GeoTouchSensor
Definition: Structs.h:4955
Multi_String
Definition: Structs.h:2468
pComponent_Geospatial
Definition: Component_Geospatial.c:307
point_XYZ
Definition: Structs.h:34
X3D_GeoViewpoint
Definition: Structs.h:5033
X3D_GeoElevationGrid
Definition: Structs.h:4719
Multi_Vec3d
Definition: Structs.h:2474
X3D_GeoMetadata
Definition: Structs.h:4836
Multi_Int32
Definition: Structs.h:2458
Uni_String
Definition: Structs.h:51
X3D_GeoTransform
Definition: Structs.h:4990
X3D_GeoLocation
Definition: Structs.h:4802
key
Definition: Viewer.h:174