35#include <libFreeWRL.h>
37#include "../vrml_parser/Structs.h"
38#include "../vrml_parser/CRoutes.h"
39#include "../main/headers.h"
41#include "../world_script/fieldSet.h"
42#include "../x3d_parser/Bindable.h"
44#include "quaternion.h"
46#include "../opengl/Frustum.h"
47#include "../opengl/Material.h"
48#include "../opengl/OpenGL_Utils.h"
49#include "../input/EAIHelpers.h"
52#include "LinearAlgebra.h"
53#include "Component_Shape.h"
54#include "Component_Geospatial.h"
56#include "../scenegraph/RenderFuncs.h"
57#include "../ui/common.h"
72void push_planetId(
int planetId);
73int current_planetId();
80static int USE_GC_FOR_LCS = TRUE;
305 int xtm_northing_first;
306 int utm_northern_hemisphere;
307 int gd_latitude_first;
314int isNodetypeGeospatial(
int nodetype,
int specversion){
317 nodetype == NODE_GeoCoordinate || nodetype == NODE_GeoElevationGrid || nodetype == NODE_GeoLOD
318 || nodetype == NODE_GeoLocation || nodetype == NODE_GeoPositionInterpolator
319 || nodetype == NODE_GeoProximitySensor || nodetype == NODE_GeoTouchSensor
320 || nodetype == NODE_GeoTransform || nodetype == NODE_GeoViewpoint || nodetype == NODE_GeoOrigin ? TRUE : FALSE;
321 if(!iret && specversion > 320) {
323 nodetype == NODE_EspduTransform || nodetype == NODE_ReceiverPdu || nodetype == NODE_SignalPdu
324 || nodetype == NODE_TransmitterPdu ? TRUE : FALSE;
328int isNodeGeospatial(
struct X3D_Node* node){
329 return isNodetypeGeospatial(node->_nodeType, X3D_PROTO(node->_executionContext)->__specversion);
335#define MF_SF_TEMPS struct Multi_Vec3d mIN; struct Multi_Vec3d mOUT; struct Multi_Vec3d gdCoords;
336#define FREE_MF_SF_TEMPS FREE_IF_NZ(gdCoords.p); FREE_IF_NZ(mOUT.p); FREE_IF_NZ(mIN.p);
338#define COMPILE_GEOSYSTEM(me) compile_geoSystem (X3D_NODE(me), me->_nodeType, &me->geoSystem, &me->__geoSystem);
340#define RADIANS_PER_DEGREE (double)0.0174532925199432957692
341#define DEGREES_PER_RADIAN (double)57.2957795130823208768
345#define UTM_SCALE (double)0.9996
346#define UTM_FALSE_EASTING 500000.0
347#define UTM_FALSE_NORTHING 10000000.0
348#define UTM_ZONE_SIZE 6.0
350#define U3TM_SCALE (double)0.9999
351#define U3TM_FALSE_EASTING 0.0
352#define U3TM_FALSE_NORTHING 0.0
353#define U3TM_ZONE_SIZE 3.0
357#define GEOEL_AA_A (double)6377563.396
358#define GEOEL_AA_F (double)299.3249646
359#define GEOEL_AM_A (double)6377340.189
360#define GEOEL_AM_F (double)299.3249646
361#define GEOEL_AN_A (double)6378160
362#define GEOEL_AN_F (double)298.25
363#define GEOEL_BN_A (double)6377483.865
364#define GEOEL_BN_F (double)299.1528128
365#define GEOEL_BR_A (double)6377397.155
366#define GEOEL_BR_F (double)299.1528128
367#define GEOEL_CC_A (double)6378206.4
368#define GEOEL_CC_F (double)294.9786982
369#define GEOEL_CD_A (double)6378249.145
370#define GEOEL_CD_F (double)293.465
371#define GEOEL_EA_A (double)6377276.345
372#define GEOEL_EA_F (double)300.8017
373#define GEOEL_EB_A (double)6377298.556
374#define GEOEL_EB_F (double)300.8017
375#define GEOEL_EC_A (double)6377301.243
376#define GEOEL_EC_F (double)300.8017
377#define GEOEL_ED_A (double)6377295.664
378#define GEOEL_ED_F (double)300.8017
379#define GEOEL_EE_A (double)6377304.063
380#define GEOEL_EE_F (double)300.8017
381#define GEOEL_EF_A (double)6377309.613
382#define GEOEL_EF_F (double)300.8017
383#define GEOEL_FA_A (double)6378155
384#define GEOEL_FA_F (double)298.3
385#define GEOEL_HE_A (double)6378200
386#define GEOEL_HE_F (double)298.3
387#define GEOEL_HO_A (double)6378270
388#define GEOEL_HO_F (double)297
389#define GEOEL_ID_A (double)6378160
390#define GEOEL_ID_F (double)298.247
391#define GEOEL_IN_A (double)6378388
392#define GEOEL_IN_F (double)297
393#define GEOEL_KA_A (double)6378245
394#define GEOEL_KA_F (double)298.3
395#define GEOEL_RF_A (double)6378137
396#define GEOEL_RF_F (double)298.257222101
397#define GEOEL_SA_A (double)6378160
398#define GEOEL_SA_F (double)298.25
399#define GEOEL_WD_A (double)6378135
400#define GEOEL_WD_F (double)298.26
401#define GEOEL_WE_A (double)6378137
402#define GEOEL_WE_F (double)298.257223563
406#define ELLIPSOIDB(typ) \
407 case typ: *semimajor = typ##_A; *flattening = 1.0/typ##_F; break;
410int nextra_ellipsoid = 1;
412int getEllipsoidParams(
int etype,
double *semimajor,
double *flattening){
415 *semimajor = *flattening = 0.0;
418 *semimajor = extra_ellipsoid[-etype].a;
419 *flattening = extra_ellipsoid[-etype].f;
445 default: printf (
"unknown ellipsoid type: %s\n", stringGEOSPATIALType(etype));
448 if(*semimajor > 0.0) iret = 1;
454#define INITIALIZE_GEOSPATIAL(me) \
455 initializeGeospatial((struct X3D_GeoOrigin **) &me->geoOrigin);
458void CONVERT_BACK_TO_GD_OR_UTMB(
Geosys *targetGeoSystem,
struct X3D_Node *GeoOrigin,
463void calculateViewingSpeed(
void);
473void clear_planets(
Stack *planet_stack){
478 for(i=0;i<vectorSize(planet_stack);i++){
479 planet = vector_get_ptr(
struct Planet,planet_stack,i);
480 if(planet && planet->gegs) deleteVector(
struct X3D_Node *, planet->gegs);
482 deleteVector(
struct Planet,planet_stack);
493 Stack *current_planet_stack;
499void *Component_Geospatial_constructor(){
504void Component_Geospatial_init(
struct tComponent_Geospatial *t){
507 t->prv = Component_Geospatial_constructor();
515 memset(&planet,0,
sizeof(
struct Planet));
516 vecset4d(planet.autoOrient.c,0.0,0.0,1.0,0.0);
517 vecsetd(planet.autoOrigin.c,0.0,0.0,0.0);
518 planet.autoOriginSet = USE_GC_FOR_LCS;
520 p->planet_stack = newStack(
struct Planet);
521 stack_push(
struct Planet,p->planet_stack,planet);
524 p->current_planet_stack = newStack(
int);
525 stack_push(
int,p->current_planet_stack,0);
526 memset(p->gcgdpars,0,50*
sizeof(
void*));
528 memset(p->fgeopars,0,50*
sizeof(
void*));
533void Component_Geospatial_clear(
struct tComponent_Geospatial *t){
538 clear_planets(p->planet_stack);
539 p->planet_stack = NULL;
541 if(p->current_planet_stack){
542 FREE_IF_NZ(p->current_planet_stack->data);
543 FREE_IF_NZ(p->current_planet_stack);
544 p->current_planet_stack = NULL;
550void push_planetId(
int planetId){
552 Stack *current_planet_stack = (
Stack*)p->current_planet_stack;
553 stack_push(
int,current_planet_stack,planetId);
555int current_planetId(){
557 Stack *current_planet_stack = (
Stack*)p->current_planet_stack;
558 return stack_top(
int,current_planet_stack);
562 Stack *current_planet_stack = (
Stack*)p->current_planet_stack;
563 stack_pop(
int,current_planet_stack);
565struct Planet *add_planet(
int planetId){
566 struct Planet planet, *ppointer;
568 memset(&planet,0,
sizeof(
struct Planet));
569 vecset4d(planet.autoOrient.c,0.0,0.0,1.0,0.0);
570 vecsetd(planet.autoOrigin.c,0.0,0.0,0.0);
571 planet.ID = planetId;
572 planet.autoOriginSet = USE_GC_FOR_LCS;
573 stack_push(
struct Planet,p->planet_stack,planet);
574 ppointer = vector_get_ptr(
struct Planet,p->planet_stack,p->planet_stack->n-1);
577struct Planet *current_planet(){
581 planetId = stack_top(
int,p->current_planet_stack);
584 for(i=0;i<vectorSize(p->planet_stack);i++){
585 planet = vector_get_ptr(
struct Planet,p->planet_stack,i);
586 if(planetId == planet->ID){
597{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},
598{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},
599{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},
600{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},
601{-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},
602{-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},
603{-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},
604{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},
605{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},
606{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},
607{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},
608{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},
609{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},
610{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},
611{-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},
612{-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},
613{-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},
614{-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},
615{-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}
618static double geoidCorrection(
double latitudeDeg,
double longitudeDeg)
629 int il, ip, il1, ip1;
630 double dl, dp, d00, d01, d10, d11, d;
632 il = (int)(longitudeDeg/10.0) + 18 -1;
633 ip = 18 - ((int)(latitudeDeg/10.0) + 9);
636 if(il1 > 35) il1 = 0;
639 d00 = (float)geoid[ip][il];
640 d01 = (float)geoid[ip1][il];
641 d10 = (float)geoid[ip][il1];
642 d11 = (float)geoid[ip1][il1];
644 dl = longitudeDeg - (-180.0 + (float)(il)*10.0);
645 dp = latitudeDeg - (90.0 - (float)(ip)*10.0);
650 d = (1.0 - dl)*(1.0 - dp)*d00 + (dl)*(1.0 - dp)*d10 + (1.0-dl)*(dp)*d01 + (dl)*(dp)*d11;
661 int geotype, lat_first, geoid;
662 double radius, flattening;
663 geotype = geoSystem->ellipsoid;
664 lat_first = geoSystem->gd_latitude_first;
665 geoid = geoSystem->geoid_height;
666 if(getEllipsoidParams(geotype,&radius,&flattening))
670 double A2 = radius*radius;
671 double F = flattening;
672 double C = A*((double)1.0 - F);
674 double Eps2 = F*((double)2.0 - F);
675 double Eps25 = (double) 0.25 * Eps2;
690 printf (
"Gd_Gc, NOT lat first\n");
691 latitude = 1; longitude = 0;
701 printf (
"Gd_Gc, have n of %d\n",inc->n);
704 for (i=0; i<n; i++) {
706 printf (
"Gd_Gc, ining lat %lf long %lf ele %lf ",LATITUDE_IN, LONGITUDE_IN, ELEVATION_IN);
709 if(geoSystem->gd_degrees == FALSE){
711 source_lat = inc[i].c[latitude];
712 source_lon = inc[i].c[longitude];
715 source_lat = RADIANS_PER_DEGREE * inc[i].c[latitude];
716 source_lon = RADIANS_PER_DEGREE * inc[i].c[longitude];
720 printf (
"Source Latitude %lf Source Longitude %lf\n",source_lat, source_lon);
723 slat = sin(source_lat);
725 clat = cos(source_lat);
728 printf (
"slat %lf slat2 %lf clat %lf\n",slat, slat2, clat);
733 Rn = A / ( (.25 - Eps25 * slat2 + .9999944354799/4) + (.25-Eps25 * slat2)/(.25 - Eps25 * slat2 + .9999944354799/4));
735 RnPh = Rn + inc[i].c[elevation];
739 double dlatin, dlongin;
740 dlatin = inc[i].c[latitude];
741 dlongin = inc[i].c[longitude];
742 if(geoSystem->gd_degrees == FALSE){
743 dlatin *= DEGREES_PER_RADIAN;
744 dlongin *= DEGREES_PER_RADIAN;
746 RnPh += geoidCorrection(dlatin,dlongin);
749 printf (
"Rn %lf RnPh %lf\n",Rn, RnPh);
752 outc[i].c[0] = RnPh * clat * cos(source_lon);
753 outc[i].c[1] = RnPh * clat * sin(source_lon);
754 outc[i].c[2] = ((C2 / A2) * Rn + inc[i].c[elevation]) * slat;
757 printf (
"Gd_Gc, outing x %lf y %lf z %lf\n", GC_X_OUT, GC_Y_OUT, GC_Z_OUT);
763double dclamp(
double fval,
double fstart,
double fend);
764static void* fwgeo_gc[50];
767 double semimajor, flattening, gd[3], gc[3];
769 getEllipsoidParams(geoSystem->ellipsoid,&semimajor,&flattening);
770 if(FALSE && flattening == 0.0){
774 veccopyd(gd,inc[i].c);
775 if(!geoSystem->gd_latitude_first) vecswizzle2d(gd);
776 if(geoSystem->gd_degrees) vecscale2d(gd,gd,RADIANS_PER_DEGREE);
777 radius = semimajor + gd[2];
778 gc[0] = cos(gd[0])*cos(gd[1])*radius;
779 gc[1] = cos(gd[0])*sin(gd[1])*radius;
780 gc[2] = sin(gd[0])*radius;
781 veccopyd(outc[i].c,gc);
787 geotype = geoSystem->ellipsoid;
788 if(geotype < 0) geotype = -geotype + GEOELLIPSOID_COUNT;
789 if(!fwgeo_gc[geotype]){
790 fwgeo_gc[geotype] = fgeo_initializeGC(semimajor,flattening);
793 veccopyd(gd,inc[i].c);
794 if(!geoSystem->gd_latitude_first) vecswizzle2d(gd);
795 if(!geoSystem->gd_degrees) vecscale2d(gd,gd,DEGREES_PER_RADIAN);
799 gd[0] = dclamp(gd[0],-90.0,90.0);
800 fgeo_gd2gc(fwgeo_gc[geotype], gd[0], gd[1], gd[2], &gc[0], &gc[1], &gc[2]);
801 veccopyd(outc[i].c,gc);
813 Gd_Gc3d_geolib(geoSystem,inc,n,outc);
821 double semimajor, flattening;
822 getEllipsoidParams(geoSystem->ellipsoid,&semimajor,&flattening);
823 if(flattening == 0.0){
826 double radius, gd[3], gc[3];
828 veccopyd(gd,inc[i].c);
829 if(!geoSystem->gd_latitude_first) vecswizzle2d(gd);
830 if(geoSystem->gd_degrees) vecscale2d(gd,gd,RADIANS_PER_DEGREE);
831 radius = semimajor + gd[2];
832 gc[0] = cos(gd[0])*cos(gd[1])*radius;
833 gc[1] = cos(gd[0])*sin(gd[1])*radius;
834 gc[2] = sin(gd[0])*radius;
835 veccopyd(outc[i].c,gc);
841 Gd_Gc3d_fw(geoSystem,inc,n,outc);
854static void Xtm_Gd3d(
Geosys *geoSystem,
struct SFVec3d *inc,
int n,
struct SFVec3d *outc,
double radius,
double flatten,
855 double scaleFactor,
double falseEasting,
double falseNorthing,
double zoneSize) {
857 int hemisphere_north,
zone, northing_first;
867 double Eccentricity = (F) * (2.0-F);
878 double longitudeOriginDeg;
879 double myeccPrimeSquared;
881 double northingDRCT1;
883 double calcConstantTerm1;
884 double calcConstantTerm2;
885 double calcConstantTerm3;
886 double calcConstantTerm4;
888 if(geoSystem->gd_latitude_first == FALSE){
893 hemisphere_north = geoSystem->utm_northern_hemisphere;
894 zone = geoSystem->xtm_zone;
895 northing_first = geoSystem->xtm_northing_first;
899 if (!northing_first) { northing = 1; easting = 0; }
904 if (!northing_first) printf (
"UTM to GD, not northing first, flipping norhting and easting\n");
908 if (northing_first) printf (
"Utm_Gd: northing first\n");
else printf (
"Utm_Gd: NOT northing_first\n");
909 if (!hemisphere_north) printf (
"Utm_Gd: NOT hemisphere_north\n");
else printf (
"Utm_Gd: hemisphere_north\n");
921 longitudeOriginDeg = (
zone -1) * zoneSize - 180. + zoneSize*.5;
922 myeccPrimeSquared = Eccentricity/(((double) 1.0) - Eccentricity);
923 eccRoot = (((double)1.0) - sqrt (((
double)1.0) - Eccentricity))/
924 (((
double)1.0) + sqrt (((
double)1.0) - Eccentricity));
926 calcConstantTerm1 = ((double)1.0) -Eccentricity/
927 ((double)4.0) - ((double)3.0) *Eccentricity*Eccentricity/
928 ((double)64.0) -((double)5.0) *Eccentricity*Eccentricity*Eccentricity/((double)256.0);
930 calcConstantTerm2 = ((double)3.0) * eccRoot/((double)2.0) - ((double)27.0) *eccRoot*eccRoot*eccRoot/((double)32.0);
931 calcConstantTerm3 = ((double)21.0) * eccRoot*eccRoot/ ((double)16.0) - ((double)55.0) *eccRoot*eccRoot*eccRoot*eccRoot/ ((double)32.0);
932 calcConstantTerm4 = ((double)151.0) *eccRoot*eccRoot*eccRoot/ ((double)96.0);
935 printf (
"zone %d\n",
zone);
936 printf (
"longitudeOriginDeg %lf\n",longitudeOriginDeg);
937 printf (
"myeccPrimeSquared %lf\n",myeccPrimeSquared);
938 printf (
"eccRoot %lf\n",eccRoot);
944 outc[i].c[elevation] = inc[i].c[elevation];
946 myEasting = inc[i].c[easting] - falseEasting;
947 if (hemisphere_north) myNorthing = inc[i].c[northing];
948 else myNorthing = inc[i].c[northing] - falseNorthing;
951 printf (
"myEasting %lf\n",myEasting);
952 printf (
"myNorthing %lf\n",myNorthing);
957 myNorthing= myNorthing / scaleFactor;
958 northingDRCT1 = myNorthing /(radius * calcConstantTerm1);
960 myphi1rad = northingDRCT1 +
961 calcConstantTerm2 * sin(((
double)2.0) *northingDRCT1)+
962 calcConstantTerm3 * sin(((
double)4.0) *northingDRCT1)+
963 calcConstantTerm4 * sin(((
double)6.0) *northingDRCT1);
965 myN1 = radius/sqrt(((
double)1.0) - Eccentricity * sin(myphi1rad) * sin (myphi1rad));
966 myT1 = tan(myphi1rad) * tan(myphi1rad);
967 myC1 = Eccentricity * cos(myphi1rad) * cos (myphi1rad);
968 myR1 = radius * (((double)1.0) - Eccentricity) / pow(((
double)1.0) - Eccentricity * sin(myphi1rad) * sin (myphi1rad), 1.5);
969 myD = myEasting/(myN1*scaleFactor);
971 Latitude = myphi1rad-(myN1*tan(myphi1rad)/myR1)*
972 (myD*myD/((
double)2.0) -
973 (((double)5.0) + ((double)3.0) *myT1+ ((double)10.0) *myC1-
974 ((double)4.0) *myC1*myC1- ((double)9.0) *myeccPrimeSquared)*
975 myD*myD*myD*myD/((
double)24.0) +(((
double)61.0) +((
double)90.0) *
976 myT1+((
double)298.0) *myC1+ ((
double)45.0) *myT1*myT1-
977 ((
double)252.0) * myeccPrimeSquared- ((
double)3.0) *myC1*myC1)*myD*myD*myD*myD*myD*myD/((double)720.0));
979 Longitude = (myD-(((double)1.0)+((double)2.0)*myT1+myC1)*myD*myD*myD/((
double)6.0)+(((
double)5.0) - ((
double)2.0) *myC1+
980 ((
double)28.0) *myT1-((
double)3.0) *myC1*myC1+
981 ((
double)8.0) *myeccPrimeSquared+((
double)24.0) *myT1*myT1)*myD*myD*myD*myD*myD/120)/cos(myphi1rad);
983 if(geoSystem->gd_degrees == FALSE){
985 outc[i].c[latitude] = Latitude ;
986 outc[i].c[longitude] = longitudeOriginDeg*RADIANS_PER_DEGREE + Longitude;
989 outc[i].c[latitude] = Latitude * DEGREES_PER_RADIAN;
990 outc[i].c[longitude] = longitudeOriginDeg + (Longitude * DEGREES_PER_RADIAN);
1006 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);
1012static void Xtm_Gd3d_geolib(
Geosys *geoSystem,
struct SFVec3d *inc,
int n,
struct SFVec3d *outc,
1013 double radius,
double flatten,
double scaleFactor,
double falseEasting,
double falseNorthing,
1016 int hemisphere_north,
zone, northing_first, geotype;
1029 double dlongitudeOrigin;
1035 if(geoSystem->gd_latitude_first == FALSE){
1040 hemisphere_north = geoSystem->utm_northern_hemisphere;
1041 zone = geoSystem->xtm_zone;
1042 northing_first = geoSystem->xtm_northing_first;
1043 geotype = geoSystem->ellipsoid;
1044 if(geotype < 0) geotype = -geotype + GEOELLIPSOID_COUNT;
1045 if(!p->fgeopars[geotype])
1046 p->fgeopars[geotype] = fgeo_initializeTM(radius, F, 1.0);
1047 fgeo = p->fgeopars[geotype];
1049 if (!northing_first) { northing = 1; easting = 0; }
1054 if (!northing_first) printf (
"UTM to GD, not northing first, flipping norhting and easting\n");
1058 if (northing_first) printf (
"Utm_Gd: northing first\n");
else printf (
"Utm_Gd: NOT northing_first\n");
1059 if (!hemisphere_north) printf (
"Utm_Gd: NOT hemisphere_north\n");
else printf (
"Utm_Gd: hemisphere_north\n");
1063 dlongitudeOrigin = (
zone -1) * zoneSize - 180. + zoneSize*.5;
1068 outc[i].c[elevation] = inc[i].c[elevation];
1070 myEasting = inc[i].c[easting] - falseEasting;
1071 if (hemisphere_north) myNorthing = inc[i].c[northing];
1072 else myNorthing = inc[i].c[northing] - falseNorthing;
1074 myEasting /= scaleFactor;
1075 myNorthing /= scaleFactor;
1078 printf (
"myEasting %lf\n",myEasting);
1079 printf (
"myNorthing %lf\n",myNorthing);
1084 fgeo_tm2gd(fgeo,myEasting, myNorthing, dlongitudeOrigin, &dLatitude, &dLongitude);
1086 if(geoSystem->gd_degrees == FALSE){
1088 outc[i].c[latitude] = dLatitude * RADIANS_PER_DEGREE ;
1089 outc[i].c[longitude] = dLongitude*RADIANS_PER_DEGREE ;
1092 outc[i].c[latitude] = dLatitude;
1093 outc[i].c[longitude] = dLongitude;
1108static void gdToUtm3d(
Geosys *geoSystem,
double *gdcoords,
double *xtmcoords);
1109static void gdTo3tm3d(
Geosys *geoSystem,
double *gdcoords,
double *xtmcoords);
1111 double semimajor, flattening;
1112 getEllipsoidParams(geoSystem->ellipsoid,&semimajor,&flattening);
1115 Xtm_Gd3d_geolib(geoSystem, inc, n, outc, semimajor, flattening, UTM_SCALE, UTM_FALSE_EASTING, UTM_FALSE_NORTHING, UTM_ZONE_SIZE);
1118 Xtm_Gd3d(geoSystem, inc, n, outc, semimajor, flattening, UTM_SCALE, UTM_FALSE_EASTING, UTM_FALSE_NORTHING, UTM_ZONE_SIZE);
1121 double xtm[3], neh[3];
1123 printf(
"in UTM_gd3d\n");
1124 printf(
"UTM y %lf x %lf h %lf zone %d\n",inc->c[0],inc->c[1],inc->c[2],geoSystem->xtm_zone);
1125 gdToUtm3d(geoSystem,outc->c, xtm);
1127 if(!geoSystem->xtm_northing_first) vecswizzle2d(neh);
1128 printf(
"UTM y %lf x %lf h %lf zone %d\n",neh[0],neh[1],neh[2],geoSystem->xtm_zone);
1130 gdTo3tm3d(geoSystem,outc->c, xtm);
1132 if(!geoSystem->xtm_northing_first) vecswizzle2d(neh);
1133 printf(
"3TM y %lf x %lf h %lf zone %d\n",neh[0],neh[1],neh[2],geoSystem->xtm_zone);
1137 double semimajor, flattening;
1138 getEllipsoidParams(geoSystem->ellipsoid,&semimajor,&flattening);
1141 Xtm_Gd3d_geolib(geoSystem, inc, n, outc, semimajor, flattening, U3TM_SCALE, U3TM_FALSE_EASTING, U3TM_FALSE_NORTHING, U3TM_ZONE_SIZE);
1144 Xtm_Gd3d(geoSystem, inc, n, outc, semimajor, flattening, U3TM_SCALE, U3TM_FALSE_EASTING, U3TM_FALSE_NORTHING, U3TM_ZONE_SIZE);
1147double getTerrainHeight(
int planetID,
Geosys *geoSystem,
struct SFVec3d *gdCoord);
1148double userHeight2ellipsoidHeight(
Geosys *geoSystem,
struct SFVec3d *gdCoord){
1149 double additionalHeight = 0.0;
1150 if(geoSystem->geoid_height == TRUE){
1154 double gddegrees[3];
1155 if(geoSystem->gd_degrees) veccopyd(gddegrees,gdCoord->c);
1156 else vecscaled(gddegrees,gdCoord->c,DEGREES_PER_RADIAN);
1157 if(!geoSystem->gd_latitude_first) vecswizzle2d(gddegrees);
1158 additionalHeight += geoidCorrection(gddegrees[0],gddegrees[1]);
1160 if(geoSystem->relativeHeight == TRUE){
1162 struct Planet *planet = current_planet();
1163 additionalHeight += getTerrainHeight( planet->ID, geoSystem, gdCoord);
1165 return additionalHeight;
1183 switch (geoSystem->spatial_system) {
1187 memcpy (gdCoords, inCoords,
sizeof (
struct SFVec3d) * n);
1188 if(geoSystem->geoid_height || geoSystem->relativeHeight)
1189 for(i=0; i < n; i++)
1190 gdCoords[i].c[2] += userHeight2ellipsoidHeight(geoSystem,&gdCoords[i]);
1193 Gd_Gc3d(geoSystem,gdCoords,n,outCoords);
1199 for (i=0; i< n; i++) {
1200 veccopyd(outCoords[i].c,inCoords[i].c);
1202 gccToGdc (geoSystem, &inCoords[i], &gdCoords[i]);
1211 Utm_Gd3d(geoSystem,inCoords,n, gdCoords);
1212 if(geoSystem->geoid_height || geoSystem->relativeHeight)
1213 for(i=0; i < n; i++)
1214 gdCoords[i].c[2] += userHeight2ellipsoidHeight(geoSystem,&gdCoords[i]);
1215 Gd_Gc3d(geoSystem,gdCoords,n,outCoords);
1223 U3tm_Gd3d(geoSystem,inCoords,n, gdCoords);
1224 if(geoSystem->geoid_height || geoSystem->relativeHeight)
1225 for(i=0; i < n; i++)
1226 gdCoords[i].c[2] += userHeight2ellipsoidHeight(geoSystem,&gdCoords[i]);
1227 Gd_Gc3d(geoSystem,gdCoords,n,outCoords);
1232 printf (
"incorrect geoSystem field, %s\n",stringGEOSPATIALType(geoSystem->spatial_system));
1240 vecdifd(outCoords[i].c,outCoords[i].c,offset->c);
1244 vrmlrot_to_quaternion(&qup,yup->c[0],yup->c[1],yup->c[2],-yup->c[3]);
1247 quaternion_rotationd(outCoords[i].c,&qup,outCoords[i].c);
1254static void initializeGeospatial (
struct X3D_GeoOrigin **nodeptr) {
1260 printf (
"\ninitializing GeoSpatial code nodeptr %u\n",*nodeptr);
1263 if (*nodeptr != NULL) {
1264 if (X3D_GEOORIGIN(*nodeptr)->_nodeType != NODE_GeoOrigin) {
1265 printf (
"expected a GeoOrigin node, but got a node of type %s\n",
1266 stringNodeType(X3D_GEOORIGIN(*nodeptr)->_nodeType));
1271 node = X3D_GEOORIGIN(*nodeptr);
1273 specversion = X3D_PROTO(node->_executionContext)->__specversion;
1276 if NODE_NEEDS_COMPILING {
1279 compile_geoSystem (X3D_NODE(node),node->_nodeType, &node->geoSystem, &node->__geoSystem);
1281 moveCoords3d(GEOSYS(node->__geoSystem), NULL, NULL,
1282 &node->geoCoords,1, &node->__movedCoords, &node->__movedgd);
1284 node->rotateYUp = FALSE;
1286 if(node->rotateYUp == TRUE)
1294 dangle = node->__movedgd.c[1];
1295 gs = GEOSYS(node->__geoSystem);
1298 dangle *= RADIANS_PER_DEGREE;
1299 dangle += RADIANS_PER_DEGREE*90.0;
1300 vrmlrot_to_quaternion (&qz,0.0, 0.0, 1.0, dangle);
1303 printf (
"GeoOrient qz angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",dangle*DEGREES_PER_RADIAN,
1304 dangle,qz.x, qz.y, qz.z,qz.w);
1307 dangle = node->__movedgd.c[0];
1308 if(gs->gd_degrees == TRUE)
1309 dangle *= RADIANS_PER_DEGREE;
1310 dangle = RADIANS_PER_DEGREE*180.0 - dangle;
1311 vrmlrot_to_quaternion (&qx,1.0, 0.0, 0.0,dangle);
1314 printf (
"GeoOrient qx angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",
1315 dangle*DEGREES_PER_RADIAN, dangle, qx.x, qx.y, qx.z,qx.w);
1319 quaternion_multiply(&qr,&qz,&qx);
1322 printf (
"GeoOrient qr %lf %lf %lf %lf\n",qr.x, qr.y, qr.z,qr.w);
1325 quaternion_to_vrmlrot(&qr, &orient.c[0], &orient.c[1], &orient.c[2], &orient.c[3]);
1327 node->__rotyup.c[i] = orient.c[i];
1331 node->__rotyup.c[i] = 0.0;
1332 node->__rotyup.c[1] = 1.0;
1336 printf (
"initializeGeospatial, __movedCoords %lf %lf %lf, ryup %d, geoSystem %d %d %d %d\n",
1337 node->__movedCoords.c[0],
1338 node->__movedCoords.c[1],
1339 node->__movedCoords.c[2],
1341 node->__geoSystem.p[0],
1342 node->__geoSystem.p[1],
1343 node->__geoSystem.p[2],
1344 node->__geoSystem.p[3]);
1345 node->__geoSystem.p[4]);
1346 node->__geoSystem.p[5]);
1347 node->__geoSystem.p[6]);
1348 node->__geoSystem.p[7]);
1349 printf (
"initializeGeospatial, done\n\n");
1360double A, F, C, A2, C2, Eps2, Eps21, Eps25, C254, C2DA, CEE,
1361 CE2, CEEps2, TwoCEE, tem, ARat1, ARat2, BRat1, BRat2, B1,B2,B3,B4,B5;
1365struct gcgd* initializeGcToGdParams(
int type,
double A,
double F) {
1368 if(type < 0) type = -type + GEOELLIPSOID_COUNT;
1369 if(p->gcgdpars[type])
return p->gcgdpars[type];
1370 g = malloc(
sizeof(
struct gcgd));
1371 p->gcgdpars[type] = g;
1376 g->C =(A) * (1-g->F);
1377 g->C2 = g->C * g->C;
1378 g->Eps2 =(g->F) * (2.0-g->F);
1379 g->Eps21 =g->Eps2 - 1.0;
1380 g->Eps25 =.25 * (g->Eps2);
1381 g->C254 =54.0 * g->C2;
1383 g->C2DA = g->C2 / A;
1384 g->CE2 = g->A2 - g->C2;
1385 g->tem = g->CE2 / g->C2;
1386 g->CEE = g->Eps2 * g->Eps2;
1387 g->TwoCEE =2.0 * g->CEE;
1388 g->CEEps2 =g->Eps2 * g->CE2;
1393 g->ARat1 =pow((A + 50005.0),2);
1394 g->ARat2 =(g->ARat1) / pow((g->C+50005.0),2);
1398 g->BRat1 =pow((A-10005.0),2);
1399 g->BRat2 =(g->BRat1) / pow((g->C-10005.0),2);
1402 g->B1=0.100225438677758E+01;
1403 g->B2=-0.393246903633930E-04;
1404 g->B3=0.241216653453483E+12;
1405 g->B4=0.133733602228679E+14;
1406 g->B5=0.984537701867943E+00;
1415 double GCC_X, GCC_Y, GCC_Z;
1417 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,
1421 printf (
"gccToGdc input %lf %lf %lf\n",GCC_X, GCC_Y, GCC_Z);
1426 if(geoSystem->gd_latitude_first == FALSE){
1427 latitude = 1; longitude = 0;
1429 getEllipsoidParams(geoSystem->ellipsoid,&A,&F);
1430 g = initializeGcToGdParams(geoSystem->ellipsoid,A,F);
1432 w2=GCC_X * GCC_X + GCC_Y * GCC_Y;
1436 testu=w2 + g->ARat2 * z2;
1437 testb=w2 + g->BRat2 * z2;
1439 if ((testb > g->BRat1) && (testu < g->ARat1))
1444 top= GCC_Z * (g->B1 + (g->B2 * w2 + g->B3) /
1445 (g->B4 + w2 * g->B5 + z2));
1461 rnn = g->A / ( (.25 - g->Eps25*s12 + .9999944354799/4) + (.25-g->Eps25*s12)/(.25 - g->Eps25*s12 + .9999944354799/4));
1469 gdc->c[elevation] = q-rnn;
1471 gdc->c[elevation] = GCC_Z / s1 + (g->Eps21 * rnn);
1472 gdc->c[latitude] = atan(top / w);
1473 gdc->c[longitude] = atan2(GCC_Y,GCC_X);
1478 wp2=GCC_X * GCC_X + GCC_Y * GCC_Y;
1482 gee=wp2 - (g->Eps21 * zp2) - g->CEEps2;
1483 alpha=cf / (gee*gee);
1484 cl=g->CEE * wp2 * alpha / gee;
1485 arg2=cl * (cl + 2.0);
1486 s1=1.0 + cl + sqrt(arg2);
1487 s=pow(s1,(1.0/3.0));
1488 p=alpha / (3.0 * pow(( s + (1.0/s) + 1.0),2));
1489 xarg= 1.0 + (g->TwoCEE * p);
1491 r2= -p * (2.0 * (1.0 - g->Eps2) * zp2 / ( q * ( 1.0 + q) ) + wp2);
1492 r1=(1.0 + (1.0 / q));
1498 ro = g->A * sqrt( .50 * (r1+r2));
1502 ro=ro - p * g->Eps2 * wp / ( 1.0 + q);
1505 arg = pow(( wp - roe),2) + zp2;
1506 v=sqrt(arg - g->Eps2 * zp2);
1507 zo=g->C2DA * GCC_Z / v;
1508 gdc->c[elevation] = sqrt(arg) * (1.0 - g->C2DA / v);
1509 top=GCC_Z+ g->tem*zo;
1510 gdc->c[latitude] = atan( top / wp );
1511 gdc->c[longitude] =atan2(GCC_Y,GCC_X);
1514 if(geoSystem->gd_degrees == TRUE){
1516 gdc->c[latitude] *= DEGREES_PER_RADIAN;
1517 gdc->c[longitude] *= DEGREES_PER_RADIAN;
1525 double gd[3],gc[3], semimajor,flattening;
1527 getEllipsoidParams(geoSystem->ellipsoid,&semimajor,&flattening);
1528 if(FALSE && flattening == 0.0){
1530 double radius, horizontal_radius;
1531 veccopyd(gc,gcc->c);
1533 radius = veclengthd(gc);
1534 horizontal_radius = veclength2d(gc);
1535 gd[0] = atan2(gc[2],horizontal_radius);
1536 gd[1] = atan2(gc[1],gc[0]);
1537 gd[2] = radius - semimajor;
1539 if(!geoSystem->gd_latitude_first) vecswizzle2d(gd);
1540 if(geoSystem->gd_degrees) vecscale2d(gd,gd,DEGREES_PER_RADIAN);
1542 veccopyd(gdc->c,gd);
1546 geotype = geoSystem->ellipsoid;
1547 if(geotype < 0) geotype = -geotype + GEOELLIPSOID_COUNT;
1548 if(!fwgeo_gc[geotype]){
1549 fwgeo_gc[geotype] = fgeo_initializeGC(semimajor,flattening);
1551 veccopyd(gc,gcc->c);
1553 fgeo_gc2gd(fwgeo_gc[geotype],gc[0],gc[1],gc[2], &gd[0],&gd[1],&gd[2]);
1554 if(!geoSystem->gd_latitude_first) vecswizzle2d(gd);
1555 if(!geoSystem->gd_degrees) vecscale2d(gd,gd,RADIANS_PER_DEGREE);
1556 veccopyd(gdc->c,gd);
1564 if(method_geolib()){
1565 gccToGdc_geolib(geoSystem,gcc,gdc);
1570 double semimajor, flattening;
1571 getEllipsoidParams(geoSystem->ellipsoid,&semimajor,&flattening);
1572 if(flattening == 0.0){
1575 double radius, horizontal_radius, gd[3], gc[3];
1576 veccopyd(gc,gcc->c);
1578 radius = veclengthd(gc);
1579 horizontal_radius = veclength2d(gc);
1580 gd[0] = atan2(gc[2],horizontal_radius);
1581 gd[1] = atan2(gc[1],gc[0]);
1582 gd[2] = radius - semimajor;
1584 if(!geoSystem->gd_latitude_first) vecswizzle2d(gd);
1585 if(geoSystem->gd_degrees) vecscale2d(gd,gd,DEGREES_PER_RADIAN);
1587 veccopyd(gdc->c,gd);
1591 gccToGdc_fw(geoSystem,gcc,gdc);
1597static void gdToXtm(
double radius,
double flattening,
double latitude,
double longitude,
double scaleFactor,
1598 double falseEasting,
double falseNorthing,
double zoneSize,
int *
zone,
double *easting,
double *northing)
1600#define DEG2RAD (PI/180.00)
1606 double longOriginradian, dlon;
1607 double eccentprime, e2, A, F;
1620 dlon = longitude * DEGREES_PER_RADIAN;
1622 *
zone = (int) (((dlon + 180.0)/zoneSize) + 1);
1624 lat_radian = latitude;
1625 long_radian = longitude;
1626 myScale = scaleFactor;
1627 longOrigin = (*
zone - 1)*zoneSize - 180.0 + zoneSize/2.0;
1628 longOriginradian = longOrigin * DEG2RAD;
1629 eccentprime = e2/(1.0-e2);
1637 NNN = A / sqrt(1.0-e2 * sin(lat_radian)*sin(lat_radian));
1638 TTT = tan(lat_radian) * tan(lat_radian);
1639 CCC = eccentprime * cos(lat_radian)*cos(lat_radian);
1640 AAA = cos(lat_radian) * (long_radian - longOriginradian);
1642 * ( ( 1. - e2/4 - 3. * e2 * e2/64
1643 - 5.0 * e2 * e2 * e2/256.0
1645 - ( 3. * e2/8 + 3. * e2 * e2/32.
1646 + 45. * e2 * e2 * e2/1024.
1647 ) * sin(2. * lat_radian)
1648 + ( 15. * e2 * e2/256. +
1649 45. * e2 * e2 * e2/1024.
1650 ) * sin(4. * lat_radian)
1651 - ( 35. * e2 * e2 * e2/3072.
1652 ) * sin(6. * lat_radian)
1657 *easting = myScale*NNN*(AAA+(1-TTT+CCC)*AAA*AAA*AAA/6
1658 + (5.-18.*TTT+TTT*TTT+72.*CCC-58.*eccentprime)*AAA*AAA*AAA*AAA*AAA/120.)
1661 *northing= myScale * ( MMM + NNN*tan(lat_radian) *
1662 ( 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.));
1664 if (latitude < 0) *northing += falseNorthing;
1667 printf (
"gdToUtm: lat %lf long %lf zone %d -> easting %lf northing %lf\n",latitude, longitude, *
zone,*easting, *northing);
1673static void gdToXtm_geolib(
int geotype,
double radius,
double flattening,
double latitude,
double longitude,
double scaleFactor,
1674 double falseEasting,
double falseNorthing,
double zoneSize,
int *
zone,
double *easting,
double *northing)
1676 double F, dlon0, dlat, dlon;
1681 if(!p->fgeopars[geotype])
1682 p->fgeopars[geotype] = fgeo_initializeTM(radius, F, 1.0);
1683 fgeo = p->fgeopars[geotype];
1687 dlon = longitude * DEGREES_PER_RADIAN;
1689 *
zone = (int) (((dlon + 180.0)/zoneSize) + 1);
1691 dlon0 = (*
zone -1) * zoneSize - 180. + zoneSize*.5;
1692 dlat = latitude * DEGREES_PER_RADIAN;
1694 fgeo_gd2tm(fgeo,dlat,dlon,dlon0,easting, northing);
1695 *easting *= scaleFactor;
1696 *northing *= scaleFactor;
1698 if (latitude < 0.0) *northing += falseNorthing;
1699 *easting += falseEasting;
1702 printf (
"gdToUtm: lat %lf long %lf zone %d -> easting %lf northing %lf\n",latitude, longitude, *
zone,*easting, *northing);
1718static void gdToUtm3d(
Geosys *geoSystem,
double *gdcoords,
double *xtmcoords) {
1719 double semimajor, flattening;
1720 double gdradians[3];
1721 int geotype, northing_first, latitude_first, is_degrees, *
zone;
1723 geotype = geoSystem->ellipsoid;
1724 northing_first = geoSystem->xtm_northing_first;
1725 latitude_first = geoSystem->gd_latitude_first;
1726 is_degrees = geoSystem->gd_degrees;
1728 if(is_degrees) vecscaled(gdradians,gdcoords,RADIANS_PER_DEGREE);
1729 else veccopyd(gdradians,gdcoords);
1730 if(!latitude_first) vecswizzle2d(gdradians);
1732 getEllipsoidParams(geotype,&semimajor,&flattening);
1733 zone = &geoSystem->xtm_zone;
1736 gdToXtm_geolib(geotype,semimajor,flattening,gdradians[0],gdradians[1], UTM_SCALE, UTM_FALSE_EASTING, UTM_FALSE_NORTHING, UTM_ZONE_SIZE,
zone, &xtmcoords[1], &xtmcoords[0]);
1739 gdToXtm(semimajor,flattening, gdradians[0],gdradians[1], UTM_SCALE, UTM_FALSE_EASTING, UTM_FALSE_NORTHING, UTM_ZONE_SIZE,
zone, &xtmcoords[1], &xtmcoords[0]);
1741 if(!northing_first) vecswizzle2d(xtmcoords);
1742 xtmcoords[2] = gdcoords[2];
1744static void gdTo3tm3d(
Geosys *geoSystem,
double *gdcoords,
double *xtmcoords) {
1745 double semimajor, flattening;
1746 double gdradians[3];
1747 int geotype, northing_first, latitude_first, is_degrees, *
zone;
1749 geotype = geoSystem->ellipsoid;
1750 northing_first = geoSystem->xtm_northing_first;
1751 latitude_first = geoSystem->gd_latitude_first;
1752 is_degrees = geoSystem->gd_degrees;
1754 if(is_degrees) vecscaled(gdradians,gdcoords,RADIANS_PER_DEGREE);
1755 else veccopyd(gdradians,gdcoords);
1756 if(!latitude_first) vecswizzle2d(gdradians);
1758 getEllipsoidParams(geotype,&semimajor,&flattening);
1759 zone = &geoSystem->xtm_zone;
1762 gdToXtm_geolib(geotype,semimajor,flattening,gdradians[0],gdradians[1], U3TM_SCALE, U3TM_FALSE_EASTING, U3TM_FALSE_NORTHING, U3TM_ZONE_SIZE,
zone, &xtmcoords[1], &xtmcoords[0]);
1765 gdToXtm(semimajor, flattening, gdradians[0],gdradians[1], U3TM_SCALE, U3TM_FALSE_EASTING, U3TM_FALSE_NORTHING, U3TM_ZONE_SIZE,
zone, &xtmcoords[1], &xtmcoords[0]);
1767 if(!northing_first) vecswizzle2d(xtmcoords);
1768 xtmcoords[2] = gdcoords[2];
1779 double dangle, gdcoords[3];
1786 if (geoSystem != NULL) {
1787 if (geoSystem->spatial_system == GEOSP_GC) {
1789 printf (
"GeoOrient - simple GC, so no orient\n");
1796 if(((
struct X3D_GeoOrigin*)geoOrigin)->rotateYUp == TRUE)
return;
1800 printf (
"GeoOrient - gdCoords->c[0,1] is %f %f\n",gdCoords->c[0],gdCoords->c[1]);
1804 veccopyd(gdcoords,gdCoords->c);
1805 if(!geoSystem->gd_latitude_first) vecswizzle2d(gdcoords);
1806 dangle = gdcoords[1];
1807 if(geoSystem->gd_degrees == TRUE)
1808 dangle *= RADIANS_PER_DEGREE;
1809 dangle += RADIANS_PER_DEGREE*90.0;
1810 vrmlrot_to_quaternion (&qz,0.0, 0.0, 1.0, dangle);
1813 printf (
"GeoOrient qz angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",((
double)90.0 + gdCoords->c[1]),
1814 RADIANS_PER_DEGREE*((
double)90.0 + gdCoords->c[1]),qz.x, qz.y, qz.z,qz.w);
1817 dangle = gdcoords[0];
1818 if(geoSystem->gd_degrees == TRUE)
1819 dangle *= RADIANS_PER_DEGREE;
1820 dangle = RADIANS_PER_DEGREE*180.0 - dangle;
1821 vrmlrot_to_quaternion (&qx,1.0, 0.0, 0.0, dangle);
1824 printf (
"GeoOrient qx angle (deg) %lf angle (rad) %lf quat: %lf %lf %lf %lf\n",
1825 ((
double)180.0 - gdCoords->c[0]), RADIANS_PER_DEGREE*((
double)180.0 - gdCoords->c[0]), qx.x, qx.y, qx.z,qx.w);
1830 quaternion_multiply(&qr, &qz, &qx);
1834 printf (
"GeoOrient qr %lf %lf %lf %lf\n",qr.x, qr.y, qr.z,qr.w);
1837 quaternion_to_vrmlrot(&qr, &orient->c[0], &orient->c[1], &orient->c[2], &orient->c[3]);
1838 vecnormald(orient->c,orient->c);
1840 printf (
"GeoOrient rotation %lf %lf %lf %lf\n",orient->c[0], orient->c[1], orient->c[2], orient->c[3]);
1858char * stringint_int2string(
struct stringint *table,
int itype){
1861 if(table[i].i == itype)
return table[i].c;
1866int stringint_string2int(
struct stringint *table,
const char *ctype){
1869 if(!strcmp(table[i].c,ctype))
return table[i].i;
1874struct stringint lookup_ellipsoids [] = {
1900struct stringint lookup_spatialreferencesys [] = {
1910 int i, specversion, nextra;
1911 indexT this_srf = INT_ID_UNDEFINED;
1912 indexT this_srf_ind = INT_ID_UNDEFINED;
1914 Geosys *srf = GEOSYS(*nodegeosys);
1917 printf (
"start of compile_geoSystem\n");
1922 srf = malloc(
sizeof(
Geosys));
1923 register_node_gc(node,(
void*)srf);
1924 *nodegeosys = X3D_NODE(srf);
1928 srf->spatial_system = GEOSP_GD;
1929 srf->ellipsoid = GEOEL_WE;
1930 srf->xtm_zone = INT_ID_UNDEFINED;
1931 srf->xtm_northing_first = TRUE;
1932 srf->utm_northern_hemisphere = TRUE;
1933 srf->gd_latitude_first = TRUE;
1934 srf->geoid_height = FALSE;
1935 specversion = X3D_PROTO(node->_executionContext)->__specversion;
1936 if(specversion > 320 && STRICT33){
1938 srf->gd_degrees = FALSE;
1941 srf->gd_degrees = TRUE;
1943 srf->relativeHeight = FALSE;
1945 if (args->n==0)
return;
1949 ee.a = ee.b = ee.f = 0.0;
1952 for (i=0; i<args->n; i++) {
1953 int itype = stringint_string2int(lookup_spatialreferencesys,args->p[i]->strptr);
1961 if (this_srf == INT_ID_UNDEFINED) {
1962 ConsoleMessage (
"geoSystem in node %s, must have GC, GD or UTM",stringNodeType(nodeType));
1966 srf->spatial_system = (int) this_srf;
1968 if (this_srf == GEOSP_GC) {
1971 }
else if (this_srf == GEOSP_GD || this_srf == GEOSP_3TM || this_srf == GEOSP_UTM) {
1972 for (i=0; i<args->n; i++) {
1973 if (i != this_srf_ind) {
1975 char *str = args->p[i]->strptr;
1977 iellipse = stringint_string2int(lookup_ellipsoids,str);
1979 srf->ellipsoid = iellipse;
1982 if (strcmp(
"latitude_first", str) == 0) {
1983 srf->gd_latitude_first = TRUE;
1984 }
else if (strcmp(
"longitude_first", str) == 0) {
1985 srf->gd_latitude_first = FALSE;
1986 }
else if(strcmp (
"WGS84",str) == 0){
1987 srf->geoid_height = TRUE;
1993 sscanf(args->p[i]->strptr,
"R%lf",&radius);
1996 }
else if (str[0] ==
'A') {
1999 sscanf(str,
"A%lf",&a);
2002 }
else if (str[0] ==
'B') {
2005 sscanf(str,
"B%lf",&b);
2008 }
else if (!strncmp(str,
"IF",2)) {
2011 sscanf(str,
"IF%lf",&invf);
2014 }
else if (str[0] ==
'F') {
2017 sscanf(str,
"F%lf",&f);
2022 if (strcmp (
"S",str) == 0) {
2023 srf->utm_northern_hemisphere = FALSE;
2024 }
else if (strcmp (
"N",str) == 0) {
2025 srf->utm_northern_hemisphere = TRUE;
2026 }
else if (str[0] ==
'Z') {
2028 sscanf(str,
"Z%d",&
zone);
2030 srf->xtm_zone =
zone;
2031 }
else if (strcmp(
"northing_first",str) == 0) {
2032 srf->xtm_northing_first = TRUE;
2033 }
else if (strcmp(
"easting_first",str) == 0) {
2034 srf->xtm_northing_first = FALSE;
2038 ConsoleMessage(
"geoSystem parameter %s not handled, in node %s",str,stringNodeType(nodeType));
2050 a = extra_ellipsoid[nextra_ellipsoid].a;
2051 b = extra_ellipsoid[nextra_ellipsoid].b;
2053 if(ee.a != 0.0 && ee.b != 0.0){
2054 ee.f = (ee.a-ee.b)/ee.a;
2055 }
else if(ee.a != 0.0){
2063 ifound = nextra_ellipsoid;
2064 for(i=1;i<nextra_ellipsoid;i++){
2065 if(extra_ellipsoid[i].a == ee.a && extra_ellipsoid[i].f == ee.f){
2070 extra_ellipsoid[ifound] = ee;
2071 srf->ellipsoid = -ifound;
2072 if(ifound == nextra_ellipsoid)
2077 printf (
"printf done compileGeoSystem\n");
2141 specversion = X3D_PROTO(node->_executionContext)->__specversion;
2143 planet = current_planet();
2144 if(!planet->autoOriginSet){
2145 if(geoOrigin && specversion < 330 ){
2148 struct SFVec3d offset, *poffset;
2153 initializeGeospatial(&geoOrigin);
2154 veccopyd(planet->autoOrigin.c,geoOrigin->__movedCoords.c);
2155 GeoOrient(X3D_NODE(geoOrigin), GEOSYS(geoOrigin->__geoSystem), &geoOrigin->__movedgd, &planet->autoOrient);
2156 planet->autoOriginSet = TRUE;
2159 user2gc(geoSystem,userCoord,1,&planet->autoOrigin);
2160 gc2gd(geoSystem,&planet->autoOrigin,1,&gdCoord);
2161 GeoOrient(X3D_NODE(geoOrigin), geoSystem, &gdCoord, &planet->autoOrient);
2162 planet->autoOriginSet = TRUE;
2170 struct SFVec3d gcCoord, lcsCoord;
2172 COMPILE_GEOSYSTEM(node)
2173 gs = GEOSYS(node->__geoSystem);
2174 FREE_IF_NZ(node->__movedCoords.p);
2175 node->__movedCoords.p = MALLOC (
struct SFVec3f *,
sizeof (
struct SFVec3f) *node->point.n);
2177 for(i=0;i<node->point.n;i++){
2178 user2gc(gs,&node->point.p[i],1,&gcCoord);
2179 gc2lcs(gs,&gcCoord,1,&lcsCoord);
2180 double2float(node->__movedCoords.p[i].c,lcsCoord.c,3);
2182 node->__movedCoords.n = node->point.n;
2195int checkX3DGeoElevationGridFields (
struct X3D_GeoElevationGrid *node,
float **points,
int *npoints) {
2209 float *texcoord = NULL;
2214 nx = node->xDimension;
2215 xSp = node->xSpacing;
2216 nz = node->zDimension;
2217 zSp = node->zSpacing;
2218 height = node->height.p;
2219 nh = node->height.n;
2221 COMPILE_GEOSYSTEM(node)
2223 if (node->__geoSystem != NULL) {
2224 mySRF = GEOSYS(node->__geoSystem)->spatial_system;
2230 rep = node->_intern;
2233 ntri = (nx && nz ? 2 * (nx-1) * (nz-1) : 0);
2239 printf (
"GeoElevationgrid: warning: x,y vs. height: %d * %d ne %d:\n", nx,nz,nh);
2241 printf (
"GeoElevationgrid: error: x,y vs. height: %d * %d ne %d:\n", nx,nz,nh);
2247 if ((nx < 2) || (nz < 2)) {
2248 printf (
"GeoElevationGrid: xDimension and zDimension less than 2 %d %d\n", nx,nz);
2256 if (!(node->texCoord)) {
2258 FREE_IF_NZ(rep->GeneratedTexCoords[0]);
2261 texcoord = rep->GeneratedTexCoords[0] = MALLOC (
float *,
sizeof (
float) * nquads * 12);
2268 newpoints = MALLOC (
float *,
sizeof (
float) * nz * nx * 3);
2270 FREE_IF_NZ(rep->actualCoord);
2271 rep->actualCoord = (
float *)newpoints;
2274 if (node->_coordIndex.n > 0) {FREE_IF_NZ(node->_coordIndex.p);}
2275 node->_coordIndex.p = MALLOC (
int *,
sizeof(
int) * nquads * 5);
2276 cindexptr = node->_coordIndex.p;
2278 node->_coordIndex.n = nquads * 5;
2280 *points = newpoints;
2281 *npoints = node->_coordIndex.n;
2284 printf (
"coordindex:\n");
2290 for (j = 0; j < (nz -1); j++) {
2291 for (i=0; i < (nx-1) ; i++) {
2293 printf (
" %d %d %d %d %d\n", j*nx+i, j*nx+i+nx, j*nx+i+nx+1, j*nx+i+1, -1);
2296#ifdef WINDING_ELEVATIONGRID
2297 *cindexptr = j*nx+i; cindexptr++;
2298 *cindexptr = j*nx+i+nx; cindexptr++;
2299 *cindexptr = j*nx+i+nx+1; cindexptr++;
2300 *cindexptr = j*nx+i+1; cindexptr++;
2301 *cindexptr = -1; cindexptr++;
2303 *cindexptr = j*nx+i; cindexptr++;
2304 *cindexptr = j*nx+i+1; cindexptr++;
2305 *cindexptr = j*nx+i+nx+1; cindexptr++;
2306 *cindexptr = j*nx+i+nx; cindexptr++;
2307 *cindexptr = -1; cindexptr++;
2315 if (!(node->texCoord)) {
2317 for (j = 0; j < (nz -1); j++) {
2318 for (i=0; i < (nx-1) ; i++) {
2320#ifdef WINDING_ELEVATIONGRID
2322 *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
2323 *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
2325 *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
2326 *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
2328 *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
2329 *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
2332 *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
2333 *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
2335 *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
2336 *texcoord = ((float)(j+1)/(nz-1)); texcoord ++;
2338 *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
2339 *texcoord = ((float)(j+0)/(nz-1)); texcoord ++;
2342 *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
2343 *texcoord = ((float) (j+0)/(nz-1)); texcoord++;
2345 *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
2346 *texcoord = ((float) (j+0)/(nz-1)); texcoord++;
2348 *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
2349 *texcoord = ((float) (j+1)/(nz-1)); texcoord++;
2352 *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
2353 *texcoord = ((float) (j+0)/(nz-1)); texcoord++;
2355 *texcoord = ((float) (i+1)/(nx-1)); texcoord++;
2356 *texcoord = ((float) (j+1)/(nz-1)); texcoord++;
2358 *texcoord = ((float) (i+0)/(nx-1)); texcoord++;
2359 *texcoord = ((float) (j+1)/(nz-1)); texcoord++;
2372 mIN.p = MALLOC (
struct SFVec3d *,
sizeof (
struct SFVec3d) * mIN.n);
2374 mOUT.n=0; mOUT.p = NULL;
2375 gdCoords.n=0; gdCoords.p = NULL;
2376 struct SFVec3d lastpoint, firstpoint;
2378 for (j=0; j<nz; j++) {
2379 for (i=0; i < nx; i++) {
2382 printf (
" %lf %lf %lf # (hei ind %d) point [%d, %d]\n",
2384 height[i+(j*nx)] * ((
double)node->yScale),
2392 if ((mySRF == GEOSP_GD) || (mySRF == GEOSP_UTM) || (mySRF == GEOSP_3TM)) {
2396 mIN.p[k].c[0] = zSp * j + node->geoGridOrigin.c[0];
2399 mIN.p[k].c[1] =xSp * i + node->geoGridOrigin.c[1];
2402 mIN.p[k].c[2] = (height[k] *(node->yScale)) + node->geoGridOrigin.c[2];
2403 veccopyd(lastpoint.c,mIN.p[k].c);
2404 if(i==0 && j==0) veccopyd(firstpoint.c,mIN.p[k].c);
2408 mIN.p[k].c[0] = zSp * j + node->geoGridOrigin.c[0];
2410 mIN.p[k].c[1] =xSp * i + node->geoGridOrigin.c[1];
2412 mIN.p[k].c[2] = (height[k] *(node->yScale)) + node->geoGridOrigin.c[2];
2420 vecprint3db(
"firstpoint",firstpoint.c,
"\n");
2421 vecprint3db(
" lastpoint",lastpoint.c,
"\n");
2422 vecprint3db(
"nx*nz",mIN.p[nx*nz -1].c,
"\n");
2423 printf (
"points before moving origin, lat, lon, height, index:\n");
2424 for (j=0; j<nz; j++) {
2425 for (i=0; i < nx; i++) {
2427 printf (
" %lf %lf %lf %d\n",mIN.p[k].c[0],
2428 mIN.p[k].c[1],mIN.p[k].c[2],k);
2438 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
2439 gs = GEOSYS(node->__geoSystem);
2441 update_origin(gs, X3D_NODE(node), &node->geoGridOrigin, X3D_GEOORIGIN(node->geoOrigin));
2444 user2gc(gs,mIN.p,mIN.n,mOUT.p);
2446 printf (
"points in gc XYZ, index:\n");
2447 for (j=0; j<nz; j++) {
2448 for (i=0; i < nx; i++) {
2449 double ci[3], co[3];
2451 veccopyd(co,mOUT.p[k].c);
2452 printf (
" %lf %lf %lf %d\n",co[0],co[1],co[2],k);
2453 veccopyd(ci,mIN.p[k].c);
2454 printf (
" %lf %lf %lf %d\n",ci[0],ci[1],ci[2],k);
2461 gc2lcs(gs,mOUT.p,mOUT.n,mOUT.p);
2468 for (j=0; j<nz; j++) {
2469 for (i=0; i < nx; i++) {
2472 double2float(newpoints,mOUT.p[k].c,3);
2477 printf (
"points converted to mesh coords, xyz index:\n");
2478 newpoints = rep->actualCoord;
2479 for (j=0; j<nz; j++) {
2480 for (i=0; i < nx; i++) {
2483 printf (
" %f %f %f %d\n",newpoints[0],newpoints[1],newpoints[2],k);
2498int planetInPlanets(
int planet,
struct Multi_Int32 *planets){
2500 for(i=0;i<planets->n;i++)
2501 if(planets->p[i] == planet) ifound = i;
2504void RegisterGeoElevationGrid(
struct X3D_Node *node,
int planetID);
2515 initializeGeospatial((
struct X3D_GeoOrigin **) &node->geoOrigin);
2516 planetID = current_planetId();
2517 COMPILE_POLY_IF_REQUIRED (NULL, NULL, node->color, node->normal, node->texCoord)
2518 CULL_FACE(node->solid)
2519 render_polyrep(node);
2520 if(!planetInPlanets(planetID,&node->__planets)){
2523 RegisterGeoElevationGrid(X3D_NODE(node), planetID);
2524 node->__planets.p = realloc(node->__planets.p,(node->__planets.n+1)*
sizeof(
int));
2525 node->__planets.p[node->__planets.n] = planetID;
2526 node->__planets.n++;
2537 struct SFVec3d gdCoord, gcCoord;
2542 planet = current_planet();
2544 printf (
"compiling GeoLocation\n");
2547 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
2548 gs = GEOSYS(node->__geoSystem);
2549 if(node->relativeHeight) gs->relativeHeight = TRUE;
2551 update_origin(gs, X3D_NODE(node), &node->geoCoords, X3D_GEOORIGIN(node->geoOrigin));
2554 MARK_SFVEC3D_INOUT_EVENT(node->geoCoords, node->__oldgeoCoords, offsetof (
struct X3D_GeoLocation, geoCoords))
2557 MARK_MFNODE_INOUT_EVENT(node->children, node->__oldChildren, offsetof (
struct X3D_GeoLocation, children))
2559 REINITIALIZE_SORTED_NODES_FIELD(node->children,node->_sortedChildren);
2568 printf (
"compiled GeoLocation\n\n");
2594 RETURN_FROM_CHILD_IF_NOT_FOR_ME
2598 prep_sibAffectors((
struct X3D_Node*)node,&node->__sibAffectors);
2606 printf (
"GeoLocation - doing normalChildren\n");
2609 normalChildren(node->children);
2612 printf (
"GeoLocation - done normalChildren\n");
2616 fin_sibAffectors((
struct X3D_Node*)node,&node->__sibAffectors);
2634 if(!renderstate()->render_vp) {
2635 geoprep(GEOSYS(node->__geoSystem),&node->geoCoords);
2638 if(renderstate()->render_boxes) extent6f_draw(node->_extent);
2646 geofin(GEOSYS(node->__geoSystem),&node->geoCoords);
2652void add_node_to_broto_context(
struct X3D_Proto *currentContext,
struct X3D_Node *node);
2654void deleteMallocedFieldValue(
int type,
union anyVrml *fieldPtr);
2658 if (childUrl->n > 0) {
2660 if (*childNode == NULL) {
2661 *childNode = createNewX3DNode(NODE_Inline);
2662 if(node->_executionContext)
2663 add_node_to_broto_context(X3D_PROTO(node->_executionContext),X3D_NODE(*childNode));
2664 ADD_PARENT(X3D_NODE(*childNode), X3D_NODE(node));
2667 deleteMallocedFieldValue(FIELDTYPE_MFString,(
union anyVrml*)&X3D_INLINE(*childNode)->url);
2668 X3D_INLINE(*childNode)->url.p = MALLOC(
struct Uni_String **,
sizeof(
struct Uni_String)*childUrl->n);
2669 for (i=0; i<childUrl->n; i++) {
2671 X3D_INLINE(*childNode)->url.p[i] = newASCIIString(childUrl->p[i]->strptr);
2674 X3D_INLINE(*childNode)->url.n = childUrl->n;
2675 X3D_INLINE(*childNode)->load = TRUE;
2679#define UNLOAD_CHILD(childNode) \
2680 if (node->childNode != NULL) \
2681 X3D_INLINE(node->childNode)->load = FALSE;
2684static void GeoLODchildren (
struct X3D_GeoLOD *node) {
2685 int load = node->__inRange;
2688 if (((node->__childloadstatus)==0) && (load)) {
2692 printf (
"GeoLODchildren - have to LOAD_CHILD for node %u (level %d)\n",node,p->geoLodLevel);
2695 LOAD_CHILD(node,&node->__child1Node,&node->child1Url);
2696 LOAD_CHILD(node,&node->__child2Node,&node->child2Url);
2697 LOAD_CHILD(node,&node->__child3Node,&node->child3Url);
2698 LOAD_CHILD(node,&node->__child4Node,&node->child4Url);
2704 node->__childloadstatus = 1;
2710static void GeoUnLODchildren (
struct X3D_GeoLOD *node) {
2711 int load = node->__inRange;
2713 if (!(load) && ((node->__childloadstatus) != 0)) {
2716 printf (
"GeoLODloadChildren, removing children from node %u level %d\n",node,p->geoLodLevel);
2718 UNLOAD_CHILD(__child1Node)
2719 UNLOAD_CHILD(__child2Node)
2720 UNLOAD_CHILD(__child3Node)
2721 UNLOAD_CHILD(__child4Node)
2723 node->__childloadstatus = 0;
2728static void GeoLODrootUrl (
struct X3D_GeoLOD *node) {
2729 int load = node->__inRange == 0;
2732 if (((node->__rooturlloadstatus)==0) && (load)) {
2734 printf (
"GeoLODrootUrl - have to LOAD_CHILD for node %u\n",node);
2737 LOAD_CHILD(node,&node->__rootUrl, &node->rootUrl);
2740 node->__rooturlloadstatus = 1;
2745static void GeoUnLODrootUrl (
struct X3D_GeoLOD *node) {
2746 int load = node->__inRange;
2748 if (!(load) && ((node->__rooturlloadstatus) != 0)) {
2750 printf (
"GeoLODloadChildren, removing rootUrl\n");
2752 node->__childloadstatus = 0;
2758void compile_GeoLOD (
struct X3D_GeoLOD * node) {
2761 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
2762 gs = GEOSYS(node->__geoSystem);
2764 update_origin(gs, X3D_NODE(node), &node->center, X3D_GEOORIGIN(node->geoOrigin));
2765 user2gc(gs,&node->center,1,&gcCoord);
2766 gc2lcs(gs,&gcCoord,1,&node->__movedCoords);
2780 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",
2784 render_vp,render_geom,render_light,render_sensitive,render_blend,render_proximity,render_collision);
2788 if (node->__level == -1) node->__level = p->geoLodLevel;
2789 else if (node->__level != p->geoLodLevel) {
2790 printf (
"hmmm - GeoLOD %p was level %d, now %d\n",node,node->__level, p->geoLodLevel);
2794 if ( node->__inRange) {
2795 printf (
"GeoLOD %u (level %d) closer\n",node,p->geoLodLevel);
2797 printf (
"GeoLOD %u (level %d) farther away\n",node,p->geoLodLevel);
2803 if (!(node->__inRange)) {
2806 GeoUnLODchildren (node);
2808 if (node->rootNode.n != 0) {
2809 for (i=0; i<node->rootNode.n; i++) {
2811 printf (
"GeoLOD %u is rendering rootNode %u",node,node->rootNode.p[i]);
2812 if (node->rootNode.p[i]!=NULL) printf (
" (%s) ",stringNodeType(X3D_NODE(node->rootNode.p[i])->_nodeType));
2816 render_node (node->rootNode.p[i]);
2818 }
else if (node->rootUrl.n != 0) {
2821 GeoLODrootUrl (node);
2824 if (node->__rootUrl != NULL) {
2826 printf (
"GeoLOD %u is rendering rootUrl %u",node,node->__rootUrl);
2827 if (node->__rootUrl != NULL) printf (
" (%s) ", stringNodeType(X3D_NODE(node->__rootUrl)->_nodeType));
2831 render_node (node->__rootUrl);
2840 GeoLODchildren (node);
2843 GeoUnLODrootUrl (node);
2846 printf (
"rendering children at %d, they are: ",p->geoLodLevel);
2847 if (node->child1Url.n>0) printf (
" :%s: ",node->child1Url.p[0]->strptr);
2848 if (node->child2Url.n>0) printf (
" :%s: ",node->child2Url.p[0]->strptr);
2849 if (node->child3Url.n>0) printf (
" :%s: ",node->child3Url.p[0]->strptr);
2850 if (node->child4Url.n>0) printf (
" :%s: ",node->child4Url.p[0]->strptr);
2856 printf (
"GeoLOD %u is rendering children %u ", node, node->__child1Node);
2857 if (node->__child1Node != NULL) printf (
" (%s) ",stringNodeType(X3D_NODE(node->__child1Node)->_nodeType));
2858 printf (
" %u ", node->__child2Node);
2859 if (node->__child2Node != NULL) printf (
" (%s) ",stringNodeType(X3D_NODE(node->__child2Node)->_nodeType));
2860 printf (
" %u ", node->__child3Node);
2861 if (node->__child3Node != NULL) printf (
" (%s) ",stringNodeType(X3D_NODE(node->__child3Node)->_nodeType));
2862 printf (
" %u ", node->__child4Node);
2863 if (node->__child4Node != NULL) printf (
" (%s) ",stringNodeType(X3D_NODE(node->__child4Node)->_nodeType));
2867 if (node->__child1Node != NULL) render_node (node->__child1Node);
2868 if (node->__child2Node != NULL) render_node (node->__child2Node);
2869 if (node->__child3Node != NULL) render_node (node->__child3Node);
2870 if (node->__child4Node != NULL) render_node (node->__child4Node);
2882 printf (
"compiling GeoMetadata\n");
2895 printf (
"compiling GeoOrigin\n");
2898 ConsoleMessage (
"compiling GeoOrigin\n");
2900 COMPILE_GEOSYSTEM(node)
2904 node->__rotyup.c[i] = 0.0;
2905 node->__rotyup.c[1] = 1.0;
2911 MARK_SFVEC3D_INOUT_EVENT(node->geoCoords,node->__oldgeoCoords,offsetof (
struct X3D_GeoOrigin, geoCoords))
2923 printf (
"compiling GeoPositionInterpolator\n");
2928 struct SFVec3d gcCoord, lcsCoord;
2929 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
2930 gs = GEOSYS(node->__geoSystem);
2931 FREE_IF_NZ(node->__movedValue.p);
2932 node->__movedValue.p = MALLOC(
struct SFVec3f*,node->keyValue.n *
sizeof(
struct SFVec3f));
2933 node->__movedValue.n = node->keyValue.n;
2934 for(i=0;i<node->keyValue.n;i++){
2935 user2gc(gs,&node->keyValue.p[i],1,&gcCoord);
2936 gc2lcs(gs,&gcCoord,1,&lcsCoord);
2937 double2float(node->__movedValue.p[i].c,lcsCoord.c,3);
2955void do_GeoPositionInterpolator (
void *innode) {
2958 int kin, kvin, counter, tmp;
2963 if (!innode)
return;
2966 if (NODE_NEEDS_COMPILING) compile_GeoPositionInterpolator(node);
2967 kvin = node->__movedValue.n;
2968 kVs_lcs = node->__movedValue.p;
2969 kVs_user = node->keyValue.p;
2975 if (node->__oldKeyValuePtr.p != node->keyValue.p) {
2977 node->__oldKeyValuePtr.p= node->keyValue.p;
2979 if (node->__oldKeyPtr.p != node->key.p) {
2981 node->__oldKeyPtr.p = node->key.p;
2986 printf(
"do_GeoPos: Position/Color interp, node %u kin %d kvin %d set_fraction %f\n",
2987 node, kin, kvin, node->set_fraction);
2991 if ((kvin == 0) || (kin == 0)) {
2992 node->value_changed.c[0] = (float) 0.0;
2993 node->value_changed.c[1] = (float) 0.0;
2994 node->value_changed.c[2] = (float) 0.0;
2995 node->geovalue_changed.c[0] = 0.0;
2996 node->geovalue_changed.c[1] = 0.0;
2997 node->geovalue_changed.c[2] = 0.0;
3001 if (kin>kvin) kin=kvin;
3004 if (node->set_fraction <= ((node->key).p[0])) {
3005 veccopyd(node->geovalue_changed.c,kVs_user[0].c);
3006 veccopy3f(node->value_changed.c,kVs_lcs[0].c);
3007 }
else if (node->set_fraction >= node->key.p[kin-1]) {
3008 memcpy ((
void *)&node->geovalue_changed, (
void *)&kVs_user[kvin-1], sizeof (
struct SFVec3d));
3009 veccopyd(node->geovalue_changed.c,kVs_user[kvin-1].c);
3010 veccopy3f(node->value_changed.c,kVs_lcs[kvin-1].c);
3013 float fpart, fdif[3];
3014 double dpart, ddif[3];
3015 counter = find_key(kin,((
float)(node->set_fraction)),node->key.p);
3018 fpart = (node->set_fraction - node->key.p[counter-1]) /
3019 (node->key.p[counter] - node->key.p[counter-1]);
3020 veclerp3f(node->value_changed.c,kVs_lcs[counter-1].c,kVs_lcs[counter].c,fpart);
3024 veclerpd(node->geovalue_changed.c,kVs_user[counter-1].c,kVs_user[counter].c,dpart);
3042 printf (
"Pos/Col, new value (%f %f %f)\n",
3043 node->value_changed.c[0],node->value_changed.c[1],node->value_changed.c[2]);
3044 printf (
"geovalue_changed %lf %lf %lf\n",node->geovalue_changed.c[0], node->geovalue_changed.c[1], node->geovalue_changed.c[2]);
3056 printf (
"compiling GeoProximitySensor\n");
3060 struct SFVec3d gcCoord, lcsCoord, gdCoord;
3061 specversion = X3D_PROTO(node->_executionContext)->__specversion;
3062 if(specversion < 330){
3065 if(!(veclengthd(node->geoCenter.c) == 0.0))
3066 veccopyd(node->center.c,node->geoCenter.c);
3069 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
3070 gs = GEOSYS(node->__geoSystem);
3071 user2gc(gs,&node->center,1,&gcCoord);
3072 gc2lcs(gs,&gcCoord,1,&node->__movedCoords);
3073 gc2gd(gs,&gcCoord,1,&gdCoord);
3074 GeoOrient(node->geoOrigin, GEOSYS(node->__geoSystem), &gdCoord, &node->__localOrient);
3076 MARK_SFVEC3D_INOUT_EVENT(node->geoCenter, node->__oldGeoCenter,offsetof (
struct X3D_GeoProximitySensor, geoCenter))
3084 printf (
"compiled GeoProximitySensor\n\n");
3095 getEllipsoidParams(gs->ellipsoid, &A, &F);
3096 if(gs->spatial_system == GEOSP_GC && veclengthd(gcCoord->c) < A/2.0){
3097 vecset4d(rotate->c,0.0,1.0,0.0,0.0);
3098 veccopyd(translate->c,gcCoord->c);
3100 gc2gd(gs,gcCoord,1, &gdCoord);
3101 tcs2gc_transform(gs,&gdCoord,rotate,translate);
3110 getEllipsoidParams(gs->ellipsoid, &A, &F);
3111 if(gs->spatial_system == GEOSP_GC && veclengthd(gcCoord->c) < A/2.0){
3113 vecset4d(rotate->c,0.0,1.0,0.0,0.0);
3114 vecsetd(translate->c,-gcCoord->c[0],-gcCoord->c[1],-gcCoord->c[2]);
3116 gc2gd(gs,gcCoord,1, &gdCoord);
3117 gc2tcs_transform(gs,&gdCoord,translate,rotate);
3126 struct SFVec3d translation1, translation2, gcCoord, gdCoord;
3127 struct SFVec4d rotation1, rotation2;
3146 getEllipsoidParams(gs->ellipsoid, &A, &F);
3147 if(gs->spatial_system == GEOSP_GC && veclengthd(userCoord->c) < A/2.0){
3150 FW_GL_TRANSLATE_D(userCoord->c[0],userCoord->c[1],userCoord->c[2]);
3152 user2gc(gs,userCoord,1,&gcCoord);
3154 gc2lcs_transform(&translation1,&rotation1);
3155 FW_GL_ROTATE_RADIANS(rotation1.c[3],rotation1.c[0],rotation1.c[1],rotation1.c[2]);
3156 FW_GL_TRANSLATE_D(translation1.c[0],translation1.c[1],translation1.c[2]);
3157 gc2gd(gs,&gcCoord,1, &gdCoord);
3158 tcs2gc_transform(gs,&gdCoord,&rotation2,&translation2);
3159 FW_GL_TRANSLATE_D(translation2.c[0],translation2.c[1],translation2.c[2]);
3160 FW_GL_ROTATE_RADIANS(rotation2.c[3],rotation2.c[0],rotation2.c[1],rotation2.c[2]);
3163 user2gc(gs,userCoord,1,&gcCoord);
3164 gc2lcs_transform(&translation1,&rotation1);
3165 FW_GL_ROTATE_RADIANS(rotation1.c[3],rotation1.c[0],rotation1.c[1],rotation1.c[2]);
3166 FW_GL_TRANSLATE_D(translation1.c[0],translation1.c[1],translation1.c[2]);
3167 tcs2gcB_transform(gs,&gcCoord,&translation2,&rotation2);
3168 FW_GL_TRANSLATE_D(translation2.c[0],translation2.c[1],translation2.c[2]);
3169 FW_GL_ROTATE_RADIANS(rotation2.c[3],rotation2.c[0],rotation2.c[1],rotation2.c[2]);
3174 if(!renderstate()->render_vp) {
3175 FW_GL_PUSH_MATRIX();
3176 geoprep0(geoSystem,userCoord);
3180void geoprepT0(
Geosys *geoSystem,
struct SFVec3d *userCoord);
3183 if(!renderstate()->render_vp) {
3186 geoprepT0(geoSystem,userCoord);
3192 if(renderstate()->render_boxes) {
3194 geoprep(GEOSYS(node->__geoSystem),&node->center);
3195 extent6f_draw(node->_extent);
3196 geofin(GEOSYS(node->__geoSystem),&node->center);
3208 static struct point_XYZ yvec = {0,0.05,0};
3209 static struct point_XYZ zvec = {0,0,-0.05};
3210 static struct point_XYZ zpvec = {0,0,0.05};
3212 struct point_XYZ t_zvec, t_yvec, t_orig, t_center;
3213 GLDOUBLE modelMatrix[16];
3214 GLDOUBLE projMatrix[16];
3215 GLDOUBLE view2prox[16];
3217 if(!((node->enabled)))
return;
3220 geoprep(GEOSYS(node->__geoSystem),&node->center);
3221 FW_GL_GETDOUBLEV(GL_MODELVIEW_MATRIX, modelMatrix);
3222 geofin(GEOSYS(node->__geoSystem),&node->center);
3223 matinverseAFFINE(view2prox,modelMatrix);
3226 transform(&t_orig,&orig,view2prox);
3228 transform(&t_center,&orig, view2prox);
3244 vecscale3f(cc,node->size.c,.5);
3245 extent6f_constructor(node->_extent,-cc[0],cc[0],-cc[1],cc[1],-cc[2],cc[2]);
3248 if(((node->size).c[0]) == 0 || ((node->size).c[1]) == 0 || ((node->size).c[2]) == 0)
return;
3250 if(fabs(cx) > ((node->size).c[0])/2 ||
3251 fabs(cy) > ((node->size).c[1])/2 ||
3252 fabs(cz) > ((node->size).c[2])/2)
return;
3260 ((node->__t1).c[0]) = (float)t_center.x;
3261 ((node->__t1).c[1]) = (float)t_center.y;
3262 ((node->__t1).c[2]) = (float)t_center.z;
3268 struct SFVec3d tcsCoord, gcCoord, userCoord;
3269 matrix_to_quaternion(&quat,modelMatrix);
3270 quaternion_normalize(&quat);
3271 quaternion_to_vrmlrot(&quat,&oo[0],&oo[1],&oo[2],&oo[3]);
3274 double2float(node->__t2.c,oo,4);
3275 float2double(tcsCoord.c,node->__t1.c,3);
3278 struct X3D_Node *boundvp = getActiveLayerBoundViewpoint();
3280 if(boundvp && boundvp->_nodeType == NODE_GeoViewpoint){
3281 struct SFVec3d gcCoord, geoCoord;
3283 user2gc(GEOSYS(gvp->__geoSystem),&gvp->position,1,&gcCoord);
3285 tcs2gc(GEOSYS(node->__geoSystem),&node->center,&tcsCoord,1,&gcCoord);
3287 gc2user(GEOSYS(node->__geoSystem),&gcCoord,1,&userCoord);
3288 veccopyd(node->__t3.c,userCoord.c);
3295void do_GeoProximitySensorTick(
void *ptr) {
3301 if (node->__oldEnabled != node->enabled) {
3302 node->__oldEnabled = node->enabled;
3305 if (!node->enabled)
return;
3310 if (!node->isActive) {
3312 printf (
"PROX - initial defaults\n");
3316 node->enterTime = TickTime();
3323 if (memcmp ((
void *) &node->position_changed,(
void *) &node->__t1,
sizeof(
struct SFColor))) {
3325 printf (
"PROX - position changed!!! \n");
3328 memcpy ((
void *) &node->position_changed,
3329 (
void *) &node->__t1,
sizeof(
struct SFColor));
3333 printf (
"do_GeoProximitySensorTick, position changed; it now is %lf %lf %lf\n",node->position_changed.c[0],
3334 node->position_changed.c[1], node->position_changed.c[2]);
3335 printf (
"nearPlane is %lf\n",Viewer.nearPlane);
3339 if(!vecsamed(node->geoCoord_changed.c,node->__t3.c)){
3340 veccopyd(node->geoCoord_changed.c,node->__t3.c);
3344 if (memcmp ((
void *) &node->orientation_changed, (
void *) &node->__t2,
sizeof(
struct SFRotation))) {
3346 printf (
"PROX - orientation changed!!!\n ");
3349 memcpy ((
void *) &node->orientation_changed,
3350 (
void *) &node->__t2,
sizeof(
struct SFRotation));
3354 if (node->isActive) {
3356 printf (
"PROX - stopping\n");
3360 node->exitTime = TickTime();
3376 printf (
"compiling GeoTouchSensor\n");
3378 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
3386void do_GeoTouchSensor (
void *ptr,
int ev,
int but1,
int over) {
3395 printf (
"%lf: TS ",TickTime());
3396 if (ev==ButtonPress) printf (
"ButtonPress ");
3397 else if (ev==ButtonRelease) printf (
"ButtonRelease ");
3398 else if (ev==KeyPress) printf (
"KeyPress ");
3399 else if (ev==KeyRelease) printf (
"KeyRelease ");
3400 else if (ev==MotionNotify) printf (
"%lf MotionNotify ");
3401 else printf (
"ev %d ",ev);
3403 if (but1) printf (
"but1 TRUE ");
else printf (
"but1 FALSE ");
3404 if (over) printf (
"over TRUE ");
else printf (
"over FALSE ");
3410 if (node->__oldEnabled != node->enabled) {
3411 node->__oldEnabled = node->enabled;
3414 if (!node->enabled)
return;
3417 if ((ev == overMark) && (over != node->isOver)) {
3419 printf (
"TS %u, isOver changed %d\n",node, over);
3421 node->isOver = over;
3427 if (ev == ButtonPress) {
3431 printf (
"touchSens %u, butPress\n",node);
3434 node->touchTime = TickTime();
3437 }
else if (ev == ButtonRelease) {
3439 printf (
"touchSens %u, butRelease\n",node);
3447 memcpy ((
void *) &node->_oldhitPoint, (
void *) &tg->RenderFuncs.ray_save_posn,
sizeof(
struct SFColor));
3450 if ((APPROX(node->_oldhitPoint.c[0],node->hitPoint_changed.c[0])!= TRUE) ||
3451 (APPROX(node->_oldhitPoint.c[1],node->hitPoint_changed.c[1])!= TRUE) ||
3452 (APPROX(node->_oldhitPoint.c[2],node->hitPoint_changed.c[2])!= TRUE)) {
3455 printf (
"GeoTouchSens, hitPoint changed: %f %f %f\n",node->hitPoint_changed.c[0],
3456 node->hitPoint_changed.c[1], node->hitPoint_changed.c[2]);
3459 memcpy ((
void *) &node->hitPoint_changed, (
void *) &node->_oldhitPoint,
sizeof(
struct SFColor));
3464 node->hitGeoCoord_changed.c[0] = (double) node->hitPoint_changed.c[0];
3465 node->hitGeoCoord_changed.c[1] = (double) node->hitPoint_changed.c[1];
3466 node->hitGeoCoord_changed.c[2] = (double) node->hitPoint_changed.c[2];
3471 printf (
"\nhitGeoCoord_changed as a GCC, %lf %lf %lf\n",
3472 node->hitGeoCoord_changed.c[0],
3473 node->hitGeoCoord_changed.c[1],
3474 node->hitGeoCoord_changed.c[2]);
3477 Geosys *gs = GEOSYS(node->__geoSystem);
3478 lcs2gc(gs,&node->hitGeoCoord_changed,1,&gcCoord);
3479 gc2user(gs,&gcCoord,1,&node->hitGeoCoord_changed);
3484 normalval.x = tg->RenderFuncs.hyp_save_norm[0];
3485 normalval.y = tg->RenderFuncs.hyp_save_norm[1];
3486 normalval.z = tg->RenderFuncs.hyp_save_norm[2];
3487 normalize_vector(&normalval);
3488 node->_oldhitNormal.c[0] = (float) normalval.x;
3489 node->_oldhitNormal.c[1] = (float) normalval.y;
3490 node->_oldhitNormal.c[2] = (float) normalval.z;
3493 MARK_SFVEC3F_INOUT_EVENT(node->hitNormal_changed,node->_oldhitNormal,offsetof (
struct X3D_GeoTouchSensor, hitNormal_changed))
3501void calculateViewingSpeedB();
3505 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
3506 gs = GEOSYS(node->__geoSystem);
3508 update_origin(gs, X3D_NODE(node), &node->position, X3D_GEOORIGIN(node->geoOrigin));
3509 user2gc(gs,&node->position,1,&gcCoord);
3510 gc2gd(gs,&gcCoord,1,&node->__movedgd);
3515 MARK_SFFLOAT_INOUT_EVENT(node->fieldOfView, node->__oldFieldOfView, offsetof (
struct X3D_GeoViewpoint, fieldOfView))
3516 MARK_SFBOOL_INOUT_EVENT(node->headlight, node->__oldHeadlight, offsetof (
struct X3D_GeoViewpoint, headlight))
3517 MARK_SFBOOL_INOUT_EVENT(node->jump, node->__oldJump, offsetof (
struct X3D_GeoViewpoint, jump))
3526 printf (
"compiled GeoViewpoint\n\n");
3529struct X3D_Node *getActiveLayerBoundViewpoint();
3530void CONVERT_BACK_TO_GD_OR_UTMC(
Geosys *targetGeoSystem,
struct X3D_Node *geoorigin,
3539 struct SFVec3d tcsCoord, gdCoord;
3543 gs = GEOSYS(node->__geoSystem);
3550 pointxyz2double(tcsCoord.c,Pos);
3551 user2gc(gs,&node->position,1,&gcCoord);
3552 gc2gd(gs,&gcCoord,1,&gdCoord);
3553 tcs2gc(gs,&gdCoord,&tcsCoord,1,&gcCoord);
3555 gc2user(gs,&gcCoord,1,&node->position);
3556 gc2gd(gs,&gcCoord,1,&gdCoord);
3563 int FREEFLY = !Viewer()->collision;
3567 struct SFVec3d translate1, translate2;
3568 struct SFVec4d rotate1, rotate2, rotate3;
3570 gc2tcs_transform(gs, &node->__movedgd, &translate1, &rotate1);
3571 gc2tcs_transform(gs, &gdCoord, &translate2, &rotate2);
3572 rotate2.c[3] = -rotate2.c[3];
3573 vrmlrot4d_to_quaternion(&q1,rotate1.c);
3574 vrmlrot4d_to_quaternion(&q2,rotate2.c);
3575 quaternion_multiply(&q3,&q1,&q2);
3576 quaternion_multiply(&q4,Quat,&q3);
3577 quaternion_to_vrmlrot4d(&q4,rotate3.c);
3578 rotate3.c[3] = -rotate3.c[3];
3579 double2float(node->orientation.c,rotate3.c,4);
3584 double oo[4], pp[3];
3585 double deltagd[3], gd[3];
3586 vecdifd(deltagd,node->__movedgd.c,gdCoord.c);
3587 veccopyd(gd,node->__movedgd.c);
3591 if(!gs->gd_latitude_first){
3593 vecswizzle2d(deltagd);
3596 if(gs->gd_degrees) {
3598 vecscale2d(deltagd,deltagd,RADIANS_PER_DEGREE);
3599 vecscale2d(gd,gd,RADIANS_PER_DEGREE);
3601 double dazimuth, dlambda;
3606 dlambda = angleNormalized(deltagd[1]);
3608 dazimuth = sin(gd[0])*dlambda;
3609 vrmlrot_to_quaternion(&qaz,0.0,1.0,0.0,-dazimuth);
3610 quaternion_multiply(&qq,Quat,&qaz);
3614 quaternion_to_vrmlrot(&qq,&oo[0],&oo[1],&oo[2],&oo[3]);
3616 double2float(node->orientation.c,oo,4);
3620 veccopyd(node->__movedgd.c,gdCoord.c);
3631 struct SFVec3d gcCoord, gdCoord, tcsCoord;
3633 gs = GEOSYS(node->__geoSystem);
3635 float2double(oo,node->orientation.c,4);
3636 vrmlrot_to_quaternion(Quat,oo[0],oo[1],oo[2], -oo[3]);
3637 user2gc(gs,&node->position,1,&gcCoord);
3638 gc2gd(gs,&gcCoord,1,&gdCoord);
3639 gc2tcs(gs,&gdCoord,&gcCoord,1,&tcsCoord);
3640 double2pointxyz(Pos,tcsCoord.c);
3657 double oo[4], gd[3];
3659 planet = current_planet();
3663 moveCoords3d(GEOSYS(node->__geoSystem),&planet->autoOrigin,&planet->autoOrient,&node->position,1,&LCpos,&node->__movedgd);
3664 vecnegated(LCpos.c,LCpos.c);
3665 double2pointxyz(Pos,LCpos.c);
3671 float2double(oo,node->orientation.c,4);
3673 vrmlrot_to_quaternion(&qoo,oo[0],oo[1],oo[2], oo[3]);
3674 GeoOrient(X3D_NODE(node->geoOrigin), GEOSYS(node->__geoSystem), &node->__movedgd, &lo);
3675 vrmlrot_to_quaternion(&qlo,lo.c[0],lo.c[1],lo.c[2], -lo.c[3]);
3676 quaternion_multiply(&qgc,&qoo,&qlo);
3679 vrmlrot_to_quaternion(&qao,planet->autoOrient.c[0],planet->autoOrient.c[1],planet->autoOrient.c[2], -planet->autoOrient.c[3]);
3680 quaternion_multiply(Quat,&qgc,&qao);
3681 quaternion_normalize(Quat);
3687 double pos[3], oo[4];
3688 struct SFVec3d GCpos, gdCoord;
3693 planet = current_planet();
3697 vrmlrot_to_quaternion(&qao,planet->autoOrient.c[0],planet->autoOrient.c[1],planet->autoOrient.c[2],planet->autoOrient.c[3]);
3698 pointxyz2double(pos,Pos);
3699 vecnegated(pos,pos);
3700 quaternion_rotationd(pos,&qao,pos);
3701 vecaddd(GCpos.c,planet->autoOrigin.c,pos);
3704 CONVERT_BACK_TO_GD_OR_UTMC(GEOSYS(node->__geoSystem), node->geoOrigin, &GCpos, &node->__movedgd, &node->position);
3709 quaternion_inverse(&qaoi,&qao);
3710 quaternion_multiply(&qgc,Quat,&qao);
3714 GeoOrient(X3D_NODE(node->geoOrigin), GEOSYS(node->__geoSystem), &node->__movedgd, &lo);
3715 vrmlrot_to_quaternion(&qlo,lo.c[0],lo.c[1],lo.c[2], lo.c[3]);
3716 quaternion_multiply(&qoo,&qgc,&qlo);
3717 quaternion_normalize(&qoo);
3718 quaternion_to_vrmlrot(&qoo,&oo[0],&oo[1],&oo[2],&oo[3]);
3720 double2float(node->orientation.c,oo,4);
3729 if (!renderstate()->render_vp)
return;
3731 if((
struct X3D_Node*)node == getActiveLayerBoundViewpoint() && !node->_donethispass){
3733 node->_donethispass = 1;
3737 printf (
"prep_GeoViewpoint called\n");
3742 struct SFVec3d translation1, translation2, gcCoord, gdCoord;
3743 struct SFVec4d rotation1, rotation2;
3745 Geosys *gs = GEOSYS(node->__geoSystem);
3751 user2gc(gs,&node->position,1,&gcCoord);
3752 gc2gd(gs,&gcCoord,1, &gdCoord);
3754 float2double(oo,node->orientation.c,4);
3756 FW_GL_ROTATE_RADIANS(oo[3],oo[0],oo[1],oo[2]);
3759 gc2tcs_transform(gs,&gdCoord,&translation2,&rotation2);
3760 FW_GL_ROTATE_RADIANS(rotation2.c[3],rotation2.c[0],rotation2.c[1],rotation2.c[2]);
3761 FW_GL_TRANSLATE_D(translation2.c[0],translation2.c[1],translation2.c[2]);
3762 lcs2gc_transform(&rotation1,&translation1);
3763 FW_GL_TRANSLATE_D(translation1.c[0],translation1.c[1],translation1.c[2]);
3764 FW_GL_ROTATE_RADIANS(rotation1.c[3],rotation1.c[0],rotation1.c[1],rotation1.c[2]);
3768 planet = current_planet();
3769 node->_prepped_planet = planet->ID;
3783 FW_GL_GETINTEGERV(GL_VIEWPORT, viewPort);
3784 if(viewPort[2] > viewPort[3]) {
3786 viewer->fieldofview = node->fieldOfView/3.1415926536*180;
3788 a1 = node->fieldOfView;
3789 a1 = atan2(sin(a1),viewPort[2]/((
float)viewPort[3]) * cos(a1));
3790 viewer->fieldofview = a1/3.1415926536*180;
3792 if(
viewer->type != VIEWER_WALK){
3794 calculateViewingSpeedB();
3795 node->_resetRelativeHeight = !node->relativeHeight;
3798 printf (
"prep_GeoViewpoint, fieldOfView %f \n",node->fieldOfView);
3804void calculateViewingSpeedB() {
3807 struct X3D_Node *boundvp = vector_back(
struct X3D_Node*,getActiveBindableStacks(tg)->viewpoint);
3809 if(boundvp->_nodeType == NODE_GeoViewpoint){
3810 double height, heightmin;
3817 COMPILE_IF_REQUIRED(X3D_NODE(node));
3818 gdCoords = &node->__movedgd;
3819 height = gdCoords->c[2];
3820 heightmin = max(abs(height),50.0);
3821 viewer->speed = heightmin * node->speedFactor;
3823 static int count = 0;
3826 printf(
"height %lf speedFactor %lf speed %lf\n",height,node->speedFactor,
viewer->speed);
3832 set_naviWidthHeightStep(
3840static void calculateExamineModeDistance(
void) {
3844Viewer()->doExamineModeDistanceCalculations = TRUE;
3851 if (!(node->isBound))
return;
3853 viewer = ViewerByLayerId(node->_layerId);
3859 if(!node->_initializedOnce) {
3860 veccopyd(node->_position.c,node->position.c);
3861 veccopy4f(node->_orientation.c,node->orientation.c);
3862 node->_initializedOnce = TRUE;
3864 if(!node->retainUserOffsets){
3865 veccopyd(node->position.c,node->_position.c);
3866 veccopy4f(node->orientation.c,node->_orientation.c);
3870 if (
viewer->transitionType != VIEWER_TRANSITION_TELEPORT &&
viewer->wasBound) {
3872 viewer->vp2rnSaved = TRUE;
3877 matcopy(
viewer->slerp_viewmatrix,bstack->viewtransformmatrix);
3878 matcopy(
viewer->slerp_posorimatrix,bstack->posorimatrix);
3882 viewer->SLERPing = FALSE;
3883 viewer->startSLERPtime = TickTime();
3885 viewer->SLERPing2 = TRUE;
3886 viewer->SLERPing2justStarted = TRUE;
3890 viewer->SLERPing = FALSE;
3891 viewer->SLERPing2 = FALSE;
3897 printf (
"bind_GeoViewpoint, setting Viewer to %lf %lf %lf orient %f %f %f %f\n",node->__movedPosition.c[0],node->__movedPosition.c[1],
3898 node->__movedPosition.c[2],node->orientation.c[0],node->orientation.c[1],node->orientation.c[2],
3899 node->orientation.c[3]);
3900 printf (
" node %u fieldOfView %f\n",node,node->fieldOfView);
3903 vrmlrot_to_quaternion (&
viewer->Quat,node->__movedOrientation.c[0],
3904 node->__movedOrientation.c[1],node->__movedOrientation.c[2],node->__movedOrientation.c[3]);
3906 calculateViewingSpeedB();
3907 node->_resetRelativeHeight = !node->relativeHeight;
3909 calculateExamineModeDistance();
3910 setMenuStatusVP (node->description->strptr);
3911 fwl_set_viewer_type (VIEWER_WALK);
3912 fwl_setCollision(TRUE);
3922 printf (
"compiling GeoLocation\n");
3926 compile_geoSystem(X3D_NODE(node),node->_nodeType,&node->geoSystem,&node->__geoSystem);
3927 gs = GEOSYS(node->__geoSystem);
3928 update_origin(gs, X3D_NODE(node), &node->geoCenter, X3D_GEOORIGIN(node->geoOrigin));
3931 MARK_SFVEC3D_INOUT_EVENT(node->geoCenter, node->__oldGeoCenter,offsetof (
struct X3D_GeoTransform, geoCenter))
3932 MARK_MFNODE_INOUT_EVENT(node->children, node->__oldChildren, offsetof (
struct X3D_GeoTransform, children))
3938 node->__do_center = verify_translate ((GLfloat *)node->center.c);
3939 node->__do_trans = verify_translate ((GLfloat *)node->translation.c);
3940 node->__do_scale = verify_scale ((GLfloat *)node->scale.c);
3941 node->__do_rotation = verify_rotate ((GLfloat *)node->rotation.c);
3942 node->__do_scaleO = verify_rotate ((GLfloat *)node->scaleOrientation.c);
3944 node->__do_anything = (node->__do_center ||
3947 node->__do_rotation ||
3950 REINITIALIZE_SORTED_NODES_FIELD(node->children,node->_sortedChildren);
3975 if(!renderstate()->render_vp) {
3977 if (node->__do_anything) {
3979 FW_GL_PUSH_MATRIX();
3982 if (node->__do_trans)
3983 FW_GL_TRANSLATE_F(node->translation.c[0],node->translation.c[1],node->translation.c[2]);
3986 if (node->__do_center)
3987 FW_GL_TRANSLATE_F(node->center.c[0],node->center.c[1],node->center.c[2]);
3990 if (node->__do_rotation) {
3991 FW_GL_ROTATE_RADIANS(node->rotation.c[3], node->rotation.c[0],node->rotation.c[1],node->rotation.c[2]);
3995 if (node->__do_scaleO) {
3996 FW_GL_ROTATE_RADIANS(node->scaleOrientation.c[3], node->scaleOrientation.c[0], node->scaleOrientation.c[1],node->scaleOrientation.c[2]);
4001 if (node->__do_scale)
4002 FW_GL_SCALE_F(node->scale.c[0],node->scale.c[1],node->scale.c[2]);
4005 if (node->__do_scaleO)
4006 FW_GL_ROTATE_RADIANS(-node->scaleOrientation.c[3], node->scaleOrientation.c[0], node->scaleOrientation.c[1],node->scaleOrientation.c[2]);
4009 if (node->__do_center)
4010 FW_GL_TRANSLATE_F(-node->center.c[0],-node->center.c[1],-node->center.c[2]);
4022 if(!renderstate()->render_vp) {
4023 if (node->__do_anything) {
4028 if((node->_renderFlags & VF_Viewpoint) == VF_Viewpoint) {
4029 FW_GL_TRANSLATE_F(((node->center).c[0]),((node->center).c[1]),((node->center).c[2])
4031 FW_GL_ROTATE_RADIANS(((node->scaleOrientation).c[3]),((node->scaleOrientation).c[0]),((node->scaleOrientation).c[1]),((node->scaleOrientation).c[2])
4033 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]))
4035 FW_GL_ROTATE_RADIANS(-(((node->scaleOrientation).c[3])),((node->scaleOrientation).c[0]),((node->scaleOrientation).c[1]),((node->scaleOrientation).c[2])
4037 FW_GL_ROTATE_RADIANS(-(((node->rotation).c[3])),((node->rotation).c[0]),((node->rotation).c[1]),((node->rotation).c[2])
4039 FW_GL_TRANSLATE_F(-(((node->center).c[0])),-(((node->center).c[1])),-(((node->center).c[2]))
4041 FW_GL_TRANSLATE_F(-(((node->translation).c[0])),-(((node->translation).c[1])),-(((node->translation).c[2]))
4046void geoprepT0(
Geosys *geoSystem,
struct SFVec3d *userCoord){
4051 struct SFVec3d translation1, translation2, gcCoord, gdCoord;
4052 struct SFVec4d rotation1, rotation2;
4061 getEllipsoidParams(gs->ellipsoid, &A, &F);
4062 if(gs->spatial_system == GEOSP_GC && veclengthd(userCoord->c) < A/2.0){
4065 FW_GL_TRANSLATE_D(-userCoord->c[0],-userCoord->c[1],-userCoord->c[2]);
4067 user2gc(gs,userCoord,1,&gcCoord);
4068 gc2gd(gs,&gcCoord,1, &gdCoord);
4070 gc2tcs_transform(gs,&gdCoord,&translation2,&rotation2);
4071 FW_GL_ROTATE_RADIANS(rotation2.c[3],rotation2.c[0],rotation2.c[1],rotation2.c[2]);
4072 FW_GL_TRANSLATE_D(translation2.c[0],translation2.c[1],translation2.c[2]);
4073 lcs2gc_transform(&rotation1,&translation1);
4074 FW_GL_TRANSLATE_D(translation1.c[0],translation1.c[1],translation1.c[2]);
4075 FW_GL_ROTATE_RADIANS(rotation1.c[3],rotation1.c[0],rotation1.c[1],rotation1.c[2]);
4078 user2gc(gs,userCoord,1,&gcCoord);
4081 gc2tcsB_transform(gs,&gcCoord,&translation2,&rotation2);
4082 FW_GL_ROTATE_RADIANS(rotation2.c[3],rotation2.c[0],rotation2.c[1],rotation2.c[2]);
4083 FW_GL_TRANSLATE_D(translation2.c[0],translation2.c[1],translation2.c[2]);
4084 lcs2gc_transform(&rotation1,&translation1);
4085 FW_GL_TRANSLATE_D(translation1.c[0],translation1.c[1],translation1.c[2]);
4086 FW_GL_ROTATE_RADIANS(rotation1.c[3],rotation1.c[0],rotation1.c[1],rotation1.c[2]);
4094 if(!renderstate()->render_vp) {
4095 FW_GL_PUSH_MATRIX();
4096 geoprepT0(geoSystem,userCoord);
4101 if(!renderstate()->render_vp) {
4104 geoprep0(geoSystem,userCoord);
4114 RETURN_FROM_CHILD_IF_NOT_FOR_ME
4133 prep_sibAffectors((
struct X3D_Node*)node,&node->__sibAffectors);
4141 printf (
"transform - doing normalChildren\n");
4143 geoprepT(GEOSYS(node->__geoSystem),&node->geoCenter);
4144 normalChildren(node->children);
4145 geofinT(GEOSYS(node->__geoSystem),&node->geoCenter);
4147 printf (
"transform - done normalChildren\n");
4151 fin_sibAffectors((
struct X3D_Node*)node,&node->__sibAffectors);
4156void CONVERT_BACK_TO_GD_OR_UTMC(
Geosys *targetGeoSystem,
struct X3D_Node *geoorigin,
4173 Geosys *geoSystem = targetGeoSystem;
4176 veccopyd(thisField->c,LCSpos->c);
4179 if (geoSystem != NULL) {
4181 if (geoSystem->spatial_system != GEOSP_GC) {
4183 gccToGdc (geoSystem, thisField, gdCoords);
4184 if(geoSystem->geoid_height || geoSystem->relativeHeight)
4185 gdCoords->c[2] -= userHeight2ellipsoidHeight(geoSystem,gdCoords);
4186 veccopyd(thisField->c,gdCoords->c);
4191 if (geoSystem->spatial_system == GEOSP_UTM || geoSystem->spatial_system == GEOSP_3TM ) {
4195 if(geoSystem->spatial_system == GEOSP_UTM){
4196 gdToUtm3d(geoSystem,thisField->c, dtemp);
4197 veccopyd(thisField->c,dtemp);
4198 }
else if(geoSystem->spatial_system == GEOSP_3TM) {
4199 gdTo3tm3d(geoSystem,thisField->c, dtemp);
4200 veccopyd(thisField->c,dtemp);
4242 float centerf[3],bottomf[3];
4243 double centerd[3], bottomd[3], spined[3], vertvecd[3];
4246 double tmin[3],tmax[3];
4249 GLDOUBLE awidth, atop, abottom, astep;
4253 nodeSystem = GEOSYS(node->__geoSystem);
4256 veccopyd(xxCoord.c,gdCoord->c);
4258 if(geoSystem->gd_latitude_first != GEOSYS(node->__geoSystem)->gd_latitude_first)
4259 vecswizzle2d(xxCoord.c);
4260 if(geoSystem->gd_degrees != GEOSYS(node->__geoSystem)->gd_degrees)
4261 if(geoSystem->gd_degrees == TRUE)
4262 vecscale2d(xxCoord.c,xxCoord.c,RADIANS_PER_DEGREE);
4264 vecscale2d(xxCoord.c,xxCoord.c,DEGREES_PER_RADIAN);
4269 if(nodeSystem->spatial_system == GEOSP_UTM || nodeSystem->spatial_system == GEOSP_3TM ) {
4273 if(nodeSystem->spatial_system == GEOSP_UTM){
4274 gdToUtm3d(GEOSYS(node->__geoSystem),xxCoord.c, dtemp);
4275 veccopyd(xxCoord.c,dtemp);
4276 }
else if(nodeSystem->spatial_system == GEOSP_3TM) {
4277 gdTo3tm3d(GEOSYS(node->__geoSystem),xxCoord.c, dtemp);
4278 veccopyd(xxCoord.c,dtemp);
4280 }
else if(nodeSystem->spatial_system == GEOSP_GC) {
4287 double emin[2],emax[2];
4289 double size[2], spacing[2];
4293 idimension[0] = node->xDimension;
4294 idimension[1] = node->zDimension;
4295 spacing[0] = node->xSpacing;
4296 spacing[1] = node->zSpacing;
4299 int swizzle_xzdimension = (XTM && nodeSystem->xtm_northing_first) || (!XTM && nodeSystem->gd_latitude_first);
4300 if(swizzle_xzdimension) {
4303 int itmp = idimension[0];
4304 idimension[0] = idimension[1];
4305 idimension[1] = itmp;
4306 vecswizzle2d(spacing);
4308 size[0] = spacing[0]*(idimension[0] -1);
4309 size[1] = spacing[1]*(idimension[1] -1);
4310 emin[0] = min(node->geoGridOrigin.c[0],node->geoGridOrigin.c[0]+size[0]);
4311 emax[0] = max(node->geoGridOrigin.c[0],node->geoGridOrigin.c[0]+size[0]);
4312 emin[1] = min(node->geoGridOrigin.c[1],node->geoGridOrigin.c[1]+size[1]);
4313 emax[1] = max(node->geoGridOrigin.c[1],node->geoGridOrigin.c[1]+size[1]);
4316 inside = xxCoord.c[0] <= emax[0] && xxCoord.c[0] >= emin[0];
4317 inside &= xxCoord.c[1] <= emax[1] && xxCoord.c[1] >= emin[1];
4320 double spinelength, vcenterd[3], pp[2];
4323 double deltah, gridpointf[3];
4327 pp[0] = xxCoord.c[0] - node->geoGridOrigin.c[0];
4328 pp[1] = xxCoord.c[1] - node->geoGridOrigin.c[1];
4331 if(swizzle_xzdimension) {
4341 double hh[4],gridheight;
4342 int i0,i1,i2,i3, ix, iz;
4343 ix = (int)(pp[0]/node->xSpacing);
4344 iz = (int)(pp[1]/node->zSpacing);
4350 i0 = iz * node->xDimension + ix;
4352 i2 = i0 + node->xDimension;
4355 hh[0] = node->height.p[i0];
4356 hh[1] = node->height.p[i1];
4357 hh[2] = node->height.p[i2];
4358 hh[3] = node->height.p[i3];
4360 cx = (pp[0] - ix*node->xSpacing)/node->xSpacing;
4361 cz = (pp[1] - iz*node->zSpacing)/node->zSpacing;
4365 gridheight = hh[0]*(1.0f - cz)*(1.0f - cx)
4366 + hh[1]*(1.0f - cz)*cx
4367 + hh[2]*cz*(1.0f - cx)
4369 gridheight *= node->yScale;
4371 *gridHeight = gridheight;
4398 float centerf[3],bottomf[3];
4399 double centerd[3], bottomd[3], spined[3], vertvecd[3], gridheight;
4400 struct SFVec3d *gdCoord, xxCoord;
4403 double tmin[3],tmax[3];
4405 GLDOUBLE awidth, atop, abottom, astep;
4407 naviinfo = (
struct sNaviInfo *)tg->Bindable.naviinfo;
4410 gdCoord = &gvp->__movedgd;
4411 geoSystem = GEOSYS(node->__geoSystem);
4412 hit = geoelevationgrid_getGDHeight0(node, gdCoord, geoSystem, &gridheight);
4413 static int count =0;
4420 if(gvp->_resetRelativeHeight){
4421 naviinfo->height = gdCoord->c[2] - gridheight;
4422 gvp->_resetRelativeHeight = FALSE;
4426 double abottom = gdCoord->c[2] - naviinfo->height;
4427 double hhh = gridheight - abottom;
4429 double hhbelowfoot = hhh;
4435 if( hhh < abottom && hhh > -fi->fallHeight)
4440 fi->hfall = hhbelowfoot;
4442 if(hhbelowfoot > fi->hfall) fi->hfall = hhbelowfoot;
4449 if( hhh >= abottom )
4454 if( fi->isClimb == 0 )
4455 fi->hclimb = hhbelowfoot;
4457 fi->hclimb = DOUBLE_MAX(fi->hclimb,hhbelowfoot);
4463 double hhabovehead = hhh - head;
4465 if( hhabovehead > 0.0 )
4475 if( fi->isClimb == 0 ){
4476 fi->hclimb = hhabovehead + abottom;
4478 fi->hclimb = DOUBLE_MAX(fi->hclimb,hhabovehead + abottom);
4506 boundvp = getActiveLayerBoundViewpoint();
4507 if(boundvp && boundvp->_nodeType == NODE_GeoViewpoint){
4509 struct Planet *planet = current_planet();;
4511 if(gvp->_prepped_planet == planet->ID) compatible = TRUE;
4513 if(node->_nodeType == NODE_GeoElevationGrid && compatible ){
4519 if(0)
if(ihit == -1){
4526double getTerrainHeight(
int planetID,
Geosys *geoSystem,
struct SFVec3d *gdCoord){
4533 if(!p->planet_stack)
return highest;
4536 for(i=0;i<vectorSize(p->planet_stack);i++){
4537 planet = vector_get_ptr(
struct Planet,p->planet_stack,i);
4538 if(planet && planet->ID == planetID) {
4539 if(!planet->gegs)
return highest;
4540 for(j=0;j<vectorSize(planet->gegs);j++){
4544 if( geoelevationgrid_getGDHeight0(geg,gdCoord,geoSystem,&gridheight) == 1){
4547 if(nfound == 1) highest = gridheight;
4548 highest = max(highest,gridheight);
4557void RegisterGeoElevationGrid(
struct X3D_Node *node,
int planetID){
4564 if(node && node->_nodeType == NODE_GeoElevationGrid){
4569 if(!p->planet_stack) p->planet_stack = newStack(
struct Planet);
4572 for(i=0;i<vectorSize(p->planet_stack);i++){
4573 planet = vector_get_ptr(
struct Planet,p->planet_stack,i);
4574 if(planet->ID == planetID) {
4582 memset(&newplanet,0,
sizeof(
struct Planet));
4583 newplanet.ID = planetID;
4585 vector_pushBack(
struct Planet,p->planet_stack,newplanet);
4586 ifound = p->planet_stack->n -1;
4587 planet = vector_get_ptr(
struct Planet,p->planet_stack,ifound);
4589 if(planet->gegs == NULL) planet->gegs = newStack(
struct X3D_Node*);
4591 vector_pushBack(
struct X3D_Node*,planet->gegs,node);
4594void unRegisterGeoElevationGrid(
struct X3D_Node *node){
4597 if(node && node->_nodeType == NODE_GeoElevationGrid){
4601 if(!p->planet_stack)
return;
4602 for(i=0;i<vectorSize(p->planet_stack);i++){
4603 struct Planet *planet = vector_get_ptr(
struct Planet,p->planet_stack,i);
4605 for(j=0;j<vectorSize(planet->gegs);j++){
4606 int k = vectorSize(planet->gegs) - j -1;
4609 vector_set(
struct X3D_Node*,planet->gegs,k,NULL);
4620 struct Planet *planet = current_planet();
4622 planet = add_planet(node->planetId);
4626 REINITIALIZE_SORTED_NODES_FIELD(node->children,node->_sortedChildren);
4637 push_planetId(node->planetId);
4639 if(!renderstate()->render_vp) {
4640 double aoo[4],ao[3];
4644 planet = current_planet();
4646 FW_GL_PUSH_MATRIX();
4647 veccopyd(ao,planet->autoOrigin.c);
4648 veccopy4d(aoo,planet->autoOrient.c);
4649 FW_GL_TRANSLATE_D(ao[0], ao[1], ao[2]);
4650 FW_GL_ROTATE_RADIANS(aoo[3], aoo[0],aoo[1],aoo[2]);
4655 if(renderstate()->render_boxes) extent6f_draw(node->_extent);
4666 RETURN_FROM_CHILD_IF_NOT_FOR_ME
4669 prep_sibAffectors((
struct X3D_Node*)node,&node->__sibAffectors);
4671 normalChildren(node->children);
4673 fin_sibAffectors((
struct X3D_Node*)node,&node->__sibAffectors);
4680 if(!renderstate()->render_vp) {
4683 if ((node->_renderFlags & VF_Viewpoint) == VF_Viewpoint) {
4684 double aoo[4],ao[3];
4687 planet = current_planet();
4688 veccopyd(ao,planet->autoOrigin.c);
4689 veccopy4d(aoo,planet->autoOrient.c);
4691 FW_GL_ROTATE_RADIANS(-aoo[3], aoo[0],aoo[1],aoo[2]);
4692 FW_GL_TRANSLATE_D(-ao[0], -ao[1], -ao[2]);
4708 moveCoords3d(geoSystem,NULL,NULL,&geo[i],1,&gc[i],&gdCoord);
4716 CONVERT_BACK_TO_GD_OR_UTMC(geoSystem,NULL,&gc[i],&gdCoord,&geo[i]);
4719void gc2lcs_transform(
struct SFVec3d *translate,
struct SFVec4d *rotate){
4724 planet = current_planet();
4725 veccopyd(translate->c,planet->autoOrigin.c);
4727 vecnegated(translate->c,translate->c);
4728 veccopy4d(rotate->c,planet->autoOrient.c);
4729 rotate->c[3] = -rotate->c[3];
4739 gc2lcs_transform(&translate, &rotate);
4740 vrmlrot_to_quaternion(&qup,rotate.c[0],rotate.c[1],rotate.c[2],rotate.c[3]);
4742 vecaddd(lcs[i].c,gc[i].c,translate.c);
4743 quaternion_rotationd(lcs[i].c,&qup,lcs[i].c);
4748void lcs2gc_transform(
struct SFVec4d *rotation,
struct SFVec3d *translation){
4753 planet = current_planet();
4754 veccopy4d(rotation->c,planet->autoOrient.c);
4755 veccopyd(translation->c,planet->autoOrigin.c);
4764 lcs2gc_transform(&rotation, &translation);
4765 vrmlrot_to_quaternion(&qup,rotation.c[0],rotation.c[1],rotation.c[2],rotation.c[3]);
4767 quaternion_rotationd(gc[i].c,&qup,lcs[i].c);
4768 vecaddd(gc[i].c,gc[i].c,translation.c);
4785 GeoOrient(NULL,geoSystem,gdcenter,rotate);
4786 rotate->c[3] = -rotate->c[3];
4787 gd2gc(geoSystem,gdcenter,1,translate);
4789 vecnegated(translate->c,translate->c);
4801 GeoOrient(NULL,geoSystem,gdcenter,rotate);
4802 gd2gc(geoSystem,gdcenter,1,translate);
4815 tcs2gc_transform(geoSystem, gdcenter, &rotate, &translate);
4816 vrmlrot_to_quaternion(&qlo,rotate.c[0],rotate.c[1],rotate.c[2], -rotate.c[3]);
4818 vecdifd(pp,gc[i].c,translate.c);
4819 quaternion_rotationd(tcs[i].c,&qlo,pp);
4832 tcs2gc_transform(geoSystem, gdcenter, &rotate, &translate);
4833 vrmlrot_to_quaternion(&qlo,rotate.c[0],rotate.c[1],rotate.c[2], rotate.c[3]);
4835 quaternion_rotationd(pp,&qlo,tcs[i].c);
4836 vecaddd(gc[i].c,translate.c,pp);
4854 tcs2gcB_transform(geoSystem, gccenter, &translate, &rotate);
4855 vrmlrot_to_quaternion(&qlo,rotate.c[0],rotate.c[1],rotate.c[2], -rotate.c[3]);
4857 vecdifd(pp,gc[i].c,translate.c);
4858 quaternion_rotationd(tcs[i].c,&qlo,pp);
4871 tcs2gcB_transform(geoSystem, gccenter, &translate, &rotate);
4872 vrmlrot_to_quaternion(&qlo,rotate.c[0],rotate.c[1],rotate.c[2], rotate.c[3]);
4874 quaternion_rotationd(pp,&qlo,tcs[i].c);
4875 vecaddd(gc[i].c,translate.c,pp);
4907 Gd_Gc3d(geoSystem,&gd[i],1,&gc[i]);
4914 gccToGdc (geoSystem, &gc[i], &gd[i]);
4917void do_GeoConvert (
void *px){
4934 if(node->__geoSystem == NULL)
4935 compile_geoSystem (X3D_NODE(node),node->_nodeType, &node->geoSystem, &node->__geoSystem);
4937 if (!vecsamed(node->__oldgeoCoords.c,node->set_geoCoords.c)) {
4940 moveCoords3d(GEOSYS(node->__geoSystem),NULL,NULL,&node->set_geoCoords,1,&node->gcCoords_changed,&gdCoord);
4941 MARK_EVENT (px, offsetof (
struct X3D_GeoConvert, gcCoords_changed));
4942 veccopyd(node->__oldgeoCoords.c,node->set_geoCoords.c);
4944 if (!vecsamed(node->__oldgcCoords.c,node->set_gcCoords.c)){
4947 CONVERT_BACK_TO_GD_OR_UTMC(GEOSYS(node->__geoSystem),NULL,&node->set_gcCoords,&gdCoord,&node->geoCoords_changed);
4948 MARK_EVENT (px, offsetof (
struct X3D_GeoConvert, geoCoords_changed));
4949 veccopyd(node->__oldgcCoords.c,node->set_gcCoords.c);