////////////////////////////////////////////////////////// /// name: Test_25_geocompass.js /// status: experimental, messy - no cleanup /// author: lukas wischounig /// license: none / open /// date: 2016 /// platform: droidscript ////////////////////////////////////////////////////////// // TODO: alocte arrays in function definions below as global arrays?? that they dont be always allocate when the function is called! var measure = 0;// global function OnStart() { // Basic settings app.SetOrientation("Portrait"); dipdir_pln_measurefrequ = "Slow"; // Sensor Measure frequency: “Fastest”, “Fast”, “Medium” or “Slow”. // TODO Function create Layout lay = app.CreateLayout("Linear", "Vertikal"); txt = app.CreateText("", 1.0, 1.0, "Multiline"); // txt = app.CreateText( text, width, height, options ); lay.AddChild(txt); app.AddLayout(lay); btn = app.CreateButton( "Measure" ); btn.SetOnTouch( GET_DIPDIR ); lay.AddChild(btn); //+++++++++++++++++++++++++++++++++++++ sns_gravity = app.CreateSensor("Gravity", dipdir_pln_measurefrequ); sns_gravity.Start(); sns_magnetic = app.CreateSensor("MagneticField", dipdir_pln_measurefrequ); sns_magnetic.Start(); MEASURE_DIPDIR(); // call function to measure while (measure == 0) { geocompass_clar = MEASURE_DIPDIR(); DIPDIR = geocompass_clar[0]; DIP = geocompass_clar[1]; GRAVITY_X_ANGLE = geocompass_clar[2]; ANG_GRAV_MAG = geocompass_clar[3]; /* txt.SetText("\n\n dipdir = " + Math.round(RAD2DEG(DIPDIR)) + "\n dip = " + Math.round(RAD2DEG(DIP)) + "\n\n\n gravity_x_angle = " + Math.round(RAD2DEG(GRAVITY_X_ANGLE)) + "\n ang_grav_mag = " + ANG_GRAV_MAG); */ /* txt.SetText("\n\n dipdir = " + Math.round(RAD2DEG(DIPDIR)) + "\n dip = " + Math.round(RAD2DEG(DIP)) + "\n\n\n gravity_x_angle = " + Math.round(RAD2DEG(GRAVITY_X_ANGLE)) + "\n ang_grav_mag = " + ANG_GRAV_MAG+ "\n\n\n quadrant= " + disptext + "\n angle_east_null_zeropoint_projected = " + angle_east_null_zeropoint_projected + "\n angle_north_null_zeropoint_projected = " + angle_north_null_zeropoint_projected) ; */ txt.SetText("\n\n dipdir = " + Math.round(RAD2DEG(DIPDIR)) + "\n dip = " + Math.round(RAD2DEG(DIP)) + "\n\n\n gravity_x_angle = " + Math.round(RAD2DEG(GRAVITY_X_ANGLE))) ; } //setInterval(MEASURE_DIPDIR, 100); // repeat function every xx milliseconds } function GET_DIPDIR() { measure = 1; } function MEASURE_DIPDIR() { // Description: Get the dip & dipdir after clar - notation for the use as a self-levelling geologcal compass. self-levelling means the funcion calculates the dip angle independent of the orientation of the measuring device in the sense, that the device can be rotated 360° ON ONE PLANE and the resulting dip-angle remains the same. // Azimut of magnetic vector only , NO TRUE NORTH CORRECTION HERE!!! // Additional output: Angle between gravity (ang_grav_mag) and magnetic vector for detection of possible local anomalies in the magnetic field (e.g. electric wires, magnetic ores, tunneling machinery etc.). Horizontal angle of the x-component of the gravity vector to align the device horizontally for the measurement of lineaments (gravity_x_angle). // The used coordinate system is a right handed system suggested by WALLBRECHER E. (1986): Tektonische und gefügeanalytische Arbeitsweisen. Ferdinand Enke Verlag Stuttgart // Input: None // Output: RESULT[0] = dipdir; RESULT[1] = dip; RESULT[2] = gravity_x_angle; RESULT[3] = ang_grav_mag; // +++++++++++++++++++++++++++++++++++++++++++ Gravity Vector and Dip - NORMALIZED var vals_grav = sns_gravity.GetValues(), gravity_r = Math.sqrt(vals_grav[1] * vals_grav[1] + vals_grav[2] * vals_grav[2] + vals_grav[3] * vals_grav[3]), gravity = []; gravity[0] = (vals_grav[2] / gravity_r) * -1; //gravity_x Notation after WALLBRECHER E. (1986) [Abb.20] gravity[1] = (vals_grav[1] / gravity_r) * -1; //gravity_y Notation after WALLBRECHER E. (1986) [Abb.20] gravity[2] = vals_grav[3] / gravity_r; //gravity_z Notation after WALLBRECHER E. (1986) [Abb.20]: +x = North, +y=East, +z = downwards //TODO: NULLWERTE KALIBRIERUNG BERÜCKSICHTIGEN var dip = Math.acos(gravity[2] / MAGNITUDE(gravity)); // dip = theta (gravity) --> see spherical coordinate systems if (dip > Math.PI / 2) {// corrcect dip (dip = Math.PI - dip); } var gravity_x_angle = Math.asin(gravity[1]);// horizontal angle of device --> for levelling device for measuring lineation if (gravity[2] < 0) { // other direction -> +pi gravity_x_angle = gravity_x_angle * -1; } // +++++++++++++++++++++++++++++++++++++++++++ Magnetic Vector - NORMALIZED var vals_magn = sns_magnetic.GetValues(), magnetic_r = Math.sqrt(vals_magn[1] * vals_magn[1] + vals_magn[2] * vals_magn[2] + vals_magn[3] * vals_magn[3]), magnetic = []; magnetic[0] = vals_magn[2] / magnetic_r; //Coord x Notation Notation after WALLBRECHER E. (1986) [Abb.20] magnetic[1] = vals_magn[1] / magnetic_r; //Coord y Notation after WALLBRECHER E. (1986) [Abb.20] magnetic[2] = (vals_magn[3] / magnetic_r) * -1; //Coord z Notation after WALLBRECHER E. (1986) [Abb.20]: +x = North, +y=East, +z = downwards var east_null = CROSSPROD(gravity, magnetic), // calculate vectorproduct (crossproduct) of gravity/magnatic vectors --> is orthogonal on both --> ALWAYS HORIZONTAL! north_null = CROSSPROD(east_null, gravity), // north_null vector calculated with crossproduct --> three orthogonal vectors (gravity; gravity/magnatic (= HORIZONTAL); & this one) zeropoint = []; // project zeropoint vector (0,0,-1) into the plane which is defined by north_null & east_null (are orthogonal on each other and always horizontal) -> get azimut of the direction in that the measurig device ist orientated in the same plane as the north_null vector zeropoint[0] = 0; zeropoint[1] = 0; zeropoint[2] = 1; var zeropoint_projected = PROJECTPOINT2PLANE(zeropoint, east_null, north_null), angle_north_null_zeropoint_projected = VECANGLE(north_null, zeropoint_projected), angle_east_null_zeropoint_projected = VECANGLE(east_null, zeropoint_projected), dipdir = 0; if ((angle_east_null_zeropoint_projected <= Math.PI && angle_east_null_zeropoint_projected > Math.PI / 2) && (angle_north_null_zeropoint_projected <= Math.PI && angle_north_null_zeropoint_projected > Math.PI / 2)) { // N->E disptext = "1 - N->E"; dipdir = Math.PI - angle_north_null_zeropoint_projected; } else if ((angle_east_null_zeropoint_projected <= Math.PI && angle_east_null_zeropoint_projected > Math.PI / 2) && (angle_north_null_zeropoint_projected <= Math.PI / 2 && angle_north_null_zeropoint_projected >= 0)) { // E->S disptext = "2 - E->S"; dipdir = (Math.PI / 2) + ((Math.PI / 2) - angle_north_null_zeropoint_projected); } else if ((angle_east_null_zeropoint_projected <= Math.PI / 2 && angle_east_null_zeropoint_projected >= 0) && (angle_north_null_zeropoint_projected <= Math.PI / 2 && angle_north_null_zeropoint_projected >= 0)) { // S->W disptext = "3 - S->W"; dipdir = Math.PI + angle_north_null_zeropoint_projected; } else { // W->N disptext = "4 - W->N"; dipdir = Math.PI + angle_north_null_zeropoint_projected; } if (gravity[2] < 0) { // other direction -> +pi dipdir = INVERTAZIMUT(dipdir); } // calculate angle between garvity vector & magnetic vector var ang_grav_mag = RAD2DEG(VECANGLE(gravity, magnetic)); var RESULT = []; RESULT[0] = dipdir; RESULT[1] = dip; RESULT[2] = gravity_x_angle; RESULT[3] = ang_grav_mag; return RESULT; } function CROSSPROD(VEC1, VEC2) { // Description: calculates the cross product of two vectors // Input: vectors VEC1 & VEC2 as arrays [x,y,z]=[0,1,2] // Output: cross product as array [x,y,z]=[0,1,2] var CROSS = []; CROSS[0] = VEC1[1] * VEC2[2] - VEC1[2] * VEC2[1]; CROSS[1] = VEC1[2] * VEC2[0] - VEC1[0] * VEC2[2]; CROSS[2] = VEC1[0] * VEC2[1] - VEC1[1] * VEC2[0]; return CROSS; } function VECANGLE(VEC1, VEC2) { // Description: calculates angle between two vectors // Input: vectors VEC1 & VEC2 as arrays [x,y,z]=[0,1,2] // Output: angle in RADIANS! return Math.acos(DOTPROD(VEC1, VEC2) / (MAGNITUDE(VEC1) * MAGNITUDE(VEC2))); } function MAGNITUDE(VEC1) { // Description: calculates the magnitude of a vector // Input: vector VEC1 as array [x,y,z]=[0,1,2] // Output: Magnitude of vector return Math.sqrt(VEC1[0] * VEC1[0] + VEC1[1] * VEC1[1] + VEC1[2] * VEC1[2]); } function DOTPROD(VEC1, VEC2) { // Description: calculates the dot product (scalar product) of two vectors // Input: vectors VEC1 & VEC2 as arrays [x,y,z]=[0,1,2] // Output: dot product (scalar) return VEC1[0] * VEC2[0] + VEC1[1] * VEC2[1] + VEC1[2] * VEC2[2]; } function VECxSCALAR(VEC, SCAL) { // Description: calculates the multiplication of a vector VEC with the scalar SCAL after // Output: vector as array [x,y,z]=[0,1,2] with RESULT_x = SCALxVEC[0], RESULT_y = SCALxVEC[1], RESULT_z = SCALxVEC[2] var RESULT = []; RESULT[0] = VEC[0] * SCAL; RESULT[1] = VEC[1] * SCAL; RESULT[2] = VEC[2] * SCAL; return RESULT; } function VECADD(VEC1, VEC2) { // Description: calculates the addition of two vectors // Output: vector as array [x,y,z]=[0,1,2] with RESULT_x = VEC1_x+VEC2_x, RESULT_y = VEC1_y+VEC2_y, RESULT_z = VEC1_z+VEC2_z var RESULT = []; RESULT[0] = VEC1[0] + VEC2[0]; RESULT[1] = VEC1[1] + VEC2[1]; RESULT[2] = VEC1[2] + VEC2[2]; return RESULT; } function PROJECTPOINT2PLANE(VECx, VECu, VECv) { // Description: calculates the ORTHOGONAL projection of a vector VECx into a plane through 0,0,0 wich is defined by the 2 ORTHOGONAL vectors VECu & VECv. (VECx is orthogonal on plane defined by VECu & VECv). // calculation after Output = ((VECx*VECu)/(VECu*VECu)*VECu+(VECx*VECv)/(VECv*VECv)*VECv) (*=dot product) after https://de.wikipedia.org/wiki/Orthogonalprojektion // Output: calculated point as array [x,y,z]=[0,1,2] return VECADD(VECxSCALAR(VECu, DOTPROD(VECx, VECu) / DOTPROD(VECu, VECu)), VECxSCALAR(VECv, DOTPROD(VECx, VECv) / DOTPROD(VECv, VECv))); } function INVERTAZIMUT(AZ) { // Description: Calculates the inverse of an azimut value (+PI) and corrects it to 0 -> Azimut -> 2*PI // Input: AZIMUT as scalar in radians // Output: AZIMUT as scalar in radians return (AZ + Math.PI) % (2 * Math.PI); } function RAD2DEG(IN) { // Description: radians to degree conversion // Input: value in radians // Output: value in degrees return IN * (180 / Math.PI); }