Commit 8d3e1fea authored by Andrey Filippov's avatar Andrey Filippov

added geomag, phone log

parent 1b048f05
<?php
// Ported from TSAGeoMag
class GeoMag{
/** The input string array which contains each line of input for the
* wmm.cof input file. Added so that all data was internal, so that
* applications do not have to mess with carrying around a data file.
* In the TSAGeoMag Class, the columns in this file are as follows:
* n, m, gnm, hnm, dgnm, dhnm
*/
private $input = array(
" 2015.0 WMM-2015 12/15/2014",
" 1 0 -29438.5 0.0 10.7 0.0",
" 1 1 -1501.1 4796.2 17.9 -26.8",
" 2 0 -2445.3 0.0 -8.6 0.0",
" 2 1 3012.5 -2845.6 -3.3 -27.1",
" 2 2 1676.6 -642.0 2.4 -13.3",
" 3 0 1351.1 0.0 3.1 0.0",
" 3 1 -2352.3 -115.3 -6.2 8.4",
" 3 2 1225.6 245.0 -0.4 -0.4",
" 3 3 581.9 -538.3 -10.4 2.3",
" 4 0 907.2 0.0 -0.4 0.0",
" 4 1 813.7 283.4 0.8 -0.6",
" 4 2 120.3 -188.6 -9.2 5.3",
" 4 3 -335.0 180.9 4.0 3.0",
" 4 4 70.3 -329.5 -4.2 -5.3",
" 5 0 -232.6 0.0 -0.2 0.0",
" 5 1 360.1 47.4 0.1 0.4",
" 5 2 192.4 196.9 -1.4 1.6",
" 5 3 -141.0 -119.4 0.0 -1.1",
" 5 4 -157.4 16.1 1.3 3.3",
" 5 5 4.3 100.1 3.8 0.1",
" 6 0 69.5 0.0 -0.5 0.0",
" 6 1 67.4 -20.7 -0.2 0.0",
" 6 2 72.8 33.2 -0.6 -2.2",
" 6 3 -129.8 58.8 2.4 -0.7",
" 6 4 -29.0 -66.5 -1.1 0.1",
" 6 5 13.2 7.3 0.3 1.0",
" 6 6 -70.9 62.5 1.5 1.3",
" 7 0 81.6 0.0 0.2 0.0",
" 7 1 -76.1 -54.1 -0.2 0.7",
" 7 2 -6.8 -19.4 -0.4 0.5",
" 7 3 51.9 5.6 1.3 -0.2",
" 7 4 15.0 24.4 0.2 -0.1",
" 7 5 9.3 3.3 -0.4 -0.7",
" 7 6 -2.8 -27.5 -0.9 0.1",
" 7 7 6.7 -2.3 0.3 0.1",
" 8 0 24.0 0.0 0.0 0.0",
" 8 1 8.6 10.2 0.1 -0.3",
" 8 2 -16.9 -18.1 -0.5 0.3",
" 8 3 -3.2 13.2 0.5 0.3",
" 8 4 -20.6 -14.6 -0.2 0.6",
" 8 5 13.3 16.2 0.4 -0.1",
" 8 6 11.7 5.7 0.2 -0.2",
" 8 7 -16.0 -9.1 -0.4 0.3",
" 8 8 -2.0 2.2 0.3 0.0",
" 9 0 5.4 0.0 0.0 0.0",
" 9 1 8.8 -21.6 -0.1 -0.2",
" 9 2 3.1 10.8 -0.1 -0.1",
" 9 3 -3.1 11.7 0.4 -0.2",
" 9 4 0.6 -6.8 -0.5 0.1",
" 9 5 -13.3 -6.9 -0.2 0.1",
" 9 6 -0.1 7.8 0.1 0.0",
" 9 7 8.7 1.0 0.0 -0.2",
" 9 8 -9.1 -3.9 -0.2 0.4",
" 9 9 -10.5 8.5 -0.1 0.3",
" 10 0 -1.9 0.0 0.0 0.0",
" 10 1 -6.5 3.3 0.0 0.1",
" 10 2 0.2 -0.3 -0.1 -0.1",
" 10 3 0.6 4.6 0.3 0.0",
" 10 4 -0.6 4.4 -0.1 0.0",
" 10 5 1.7 -7.9 -0.1 -0.2",
" 10 6 -0.7 -0.6 -0.1 0.1",
" 10 7 2.1 -4.1 0.0 -0.1",
" 10 8 2.3 -2.8 -0.2 -0.2",
" 10 9 -1.8 -1.1 -0.1 0.1",
" 10 10 -3.6 -8.7 -0.2 -0.1",
" 11 0 3.1 0.0 0.0 0.0",
" 11 1 -1.5 -0.1 0.0 0.0",
" 11 2 -2.3 2.1 -0.1 0.1",
" 11 3 2.1 -0.7 0.1 0.0",
" 11 4 -0.9 -1.1 0.0 0.1",
" 11 5 0.6 0.7 0.0 0.0",
" 11 6 -0.7 -0.2 0.0 0.0",
" 11 7 0.2 -2.1 0.0 0.1",
" 11 8 1.7 -1.5 0.0 0.0",
" 11 9 -0.2 -2.5 0.0 -0.1",
" 11 10 0.4 -2.0 -0.1 0.0",
" 11 11 3.5 -2.3 -0.1 -0.1",
" 12 0 -2.0 0.0 0.1 0.0",
" 12 1 -0.3 -1.0 0.0 0.0",
" 12 2 0.4 0.5 0.0 0.0",
" 12 3 1.3 1.8 0.1 -0.1",
" 12 4 -0.9 -2.2 -0.1 0.0",
" 12 5 0.9 0.3 0.0 0.0",
" 12 6 0.1 0.7 0.1 0.0",
" 12 7 0.5 -0.1 0.0 0.0",
" 12 8 -0.4 0.3 0.0 0.0",
" 12 9 -0.4 0.2 0.0 0.0",
" 12 10 0.2 -0.9 0.0 0.0",
" 12 11 -0.9 -0.2 0.0 0.0",
" 12 12 0.0 0.7 0.0 0.0",
);
/**
* Geodetic altitude in km. An input,
* but set to zero in this class. Changed
* back to an input in version 5. If not specified,
* then is 0.
*/
private $alt = 0; // double
/**
* Geodetic latitude in deg. An input.
*/
private $glat = 0; // double
/**
* Geodetic longitude in deg. An input.
*/
private $glon = 0; // double
/**
* Time in decimal years. An input.
*/
private $time = 0; // double
/**
* Geomagnetic declination in deg.
* East is positive, West is negative.
* (The negative of variation.)
*/
private $dec = 0; // double
/**
* Geomagnetic inclination in deg.
* Down is positive, up is negative.
*/
private $dip = 0; // double
/**
* Geomagnetic total intensity, in nano Teslas.
*/
private $ti = 0; // double
/**
* Geomagnetic grid variation, referenced to
* grid North. Not calculated or output in version 5.0.
*/
//private double gv = 0;
/**
* The maximum number of degrees of the spherical harmonic model.
*/
private $maxdeg = 12; // int
/**
* The maximum order of spherical harmonic model.
*/
private $maxord; // int
/** Added in version 5. In earlier versions the date for the calculation was held as a
* constant. Now the default date is set to 2.5 years plus the epoch read from the
* input file.
*/
private $defaultDate = 2017.5; // double
/** Added in version 5. In earlier versions the altitude for the calculation was held as a
* constant at 0. In version 5, if no altitude is specified in the calculation, this
* altitude is used by default.
*/
private $defaultAltitude = 0; // double
/**
* The Gauss coefficients of main geomagnetic model (nt).
*/
// private double c[][] = new double[13][13]; // double // array_fill(0, $rows * $cols, 0.0);
private $c = array();
// for ($i = 0; $i<13; $i++) $c[]= array_fill(0, 13, 0.0);
/**
* The Gauss coefficients of secular geomagnetic model (nt/yr).
*/
private $cd = array(); // new double[13][13];
/**
* The time adjusted geomagnetic gauss coefficients (nt).
*/
private $tc = array(); // new double[13][13];
/**
* The theta derivative of p(n,m) (unnormalized).
*/
private $dp = array(); // new double[13][13];
/**
* The Schmidt normalization factors.
*/
private $snorm = array(); // new double[169];
/**
* The sine of (m*spherical coord. longitude).
*/
private $sp = array(); // new double[13];
/**
* The cosine of (m*spherical coord. longitude).
*/
private $cp = array(); // new double[13];
private $fn = array(); // new double[13];
private $fm = array(); // new double[13];
/**
* The associated Legendre polynomials for m=1 (unnormalized).
*/
private $pp = array(); // new double[13];
private $k = array(); // new double[13][13];
/**
* The variables otime (old time), oalt (old altitude),
* olat (old latitude), olon (old longitude), are used to
* store the values used from the previous calculation to
* save on calculation time if some inputs don't change.
*/
private $otime, $oalt, $olat, $olon; // double
/** The date in years, for the start of the valid time of the fit coefficients */
private $epoch; // double
/** bx is the north south field intensity
* by is the east west field intensity
* bz is the vertical field intensity positive downward
* bh is the horizontal field intensity
*/
private $bx, $by, $bz, $bh; // double
private $re, $a2, $b2, $c2, $a4, $b4, $c4; // double
private $r, $d, $ca, $sa, $ct, $st; // double // even though these only occur in one method, they must be
// created here, or won't have correct values calculated
// These values are only recalculated if the altitude changes.
public function __construct()
{
// initialize arrays [13][13]
for ($i = 0; $i<13; $i++) $this->c [] = array_fill(0, 13, 0.0);
for ($i = 0; $i<13; $i++) $this->cd[] = array_fill(0, 13, 0.0);
for ($i = 0; $i<13; $i++) $this->tc[] = array_fill(0, 13, 0.0);
for ($i = 0; $i<13; $i++) $this->dp[] = array_fill(0, 13, 0.0);
$this->snorm = array_fill(0, 169, 0.0);
$this->sp = array_fill(0, 13, 0.0);
$this->cp = array_fill(0, 13, 0.0);
$this->fn = array_fill(0, 13, 0.0);
$this->fm = array_fill(0, 13, 0.0);
$this->pp = array_fill(0, 13, 0.0);
for ($i = 0; $i<13; $i++) $this->k[] = array_fill(0, 13, 0.0);
//read model data from file and initialize the GeoMag routine
$this->initModel();
}
private function initModel()
{
$this->glat = 0;
$this->glon = 0;
// INITIALIZE CONSTANTS
$this->maxord = $this->maxdeg;
$this->sp[0] = 0.0;
$this->cp[0] = $this->snorm[0] = $this->pp[0] = 1.0;
$this->dp[0][0] = 0.0;
/**
* Semi-major axis of WGS-84 ellipsoid, in km.
*/
$a = 6378.137; // double
/**
* Semi-minor axis of WGS-84 ellipsoid, in km.
*/
$b = 6356.7523142; // double
/**
* Mean radius of IAU-66 ellipsoid, in km.
*/
$this->re = 6371.2;
$this->a2 = $a * $a;
$this->b2 = $b * $b;
$this->c2 = $this->a2 - $this->b2;
$this->a4 = $this->a2 * $this->a2;
$this->b4 = $this->b2 * $this->b2;
$this->c4 = $this->a4 - $this->b4;
// Skipping file read, using hardwired data only
$this->setCoeff();
// CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED
$this->snorm[0] = 1.0;
for ($n = 1; $n <= $this->maxord; $n++){
$this->snorm[$n] = $this->snorm[$n - 1] * (2 * $n - 1) / $n;
$j = 2;
for($m = 0, $D1 = 1, $D2 = ($n - $m + $D1) / $D1; $D2 > 0; $D2--, $m += $D1){
$this->k[$m][$n] = ((($n - 1) * ($n - 1))-($m * $m)) / ((2 * $n - 1) * (2 * $n-3));
if($m > 0){
$flnmj = (($n - $m + 1) * $j) / ($n + $m);
$this->snorm[$n + $m * 13] = $this->snorm[$n + ($m -1 ) * 13] * sqrt($flnmj);
$j = 1;
$this->c [$n][$m-1] = $this->snorm[$n + $m * 13] * $this->c [$n][$m-1];
$this->cd[$n][$m-1] = $this->snorm[$n + $m * 13] * $this->cd[$n][$m-1];
}
$this->c [$m][$n] = $this->snorm[$n + $m * 13] * $this->c [$m][$n];
$this->cd[$m][$n] = $this->snorm[$n + $m * 13] * $this->cd[$m][$n];
} //for(m...)
$this->fn[$n] = ($n+1);
$this->fm[$n] = $n;
} //for(n...)
$this->k[1][1] = 0.0;
$this->otime = $this->oalt = $this->olat = $this->olon = -1000.0;
}
/*
* @param fLat The latitude in decimal degrees.
* @param fLon The longitude in decimal degrees.
* @param year The date as a decimal year.
* @param altitude The altitude in kilometers.
*/
private function calcGeoMag($fLat, $fLon, $year=null, $altitude=null)
{
if ($year === null) $year= $this->defaultDate;
if ($altitude === null) $altitude= $this->defaultAltitude;
$this->glat = $fLat;
$this->glon = $fLon;
$this->alt = $altitude;
/**
* The date in decimal years for calculating the magnetic field components.
*/
$this->time = $year;
$dt = $this->time - $this->epoch;
$pi = pi();;
$dtr = ($pi/180.0);
$rlon = $this->glon * $dtr;
$rlat = $this->glat * $dtr;
$srlon = sin($rlon);
$srlat = sin($rlat);
$crlon = cos($rlon);
$crlat = cos($rlat);
$srlat2 = $srlat * $srlat;
$crlat2 = $crlat * $crlat;
$this->sp[1] = $srlon;
$this->cp[1] = $crlon;
// CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS.
if ($this->alt != $this->oalt || $this->glat != $this->olat){
$q = sqrt($this->a2 - $this->c2 * $srlat2);
$q1 = $this->alt * $q;
$q2 = (($q1 + $this->a2) / ($q1 + $this->b2)) * (($q1 + $this->a2) / ($q1 + $this->b2));
$this->ct = $srlat / sqrt($q2 * $crlat2 + $srlat2);
$this->st = sqrt(1.0 - ($this->ct * $this->ct));
$r2 = (($this->alt * $this->alt) + 2.0 * $q1 + ($this->a4 - $this->c4 * $srlat2) / ($q * $q));
$this->r = sqrt($r2);
$this->d = sqrt($this->a2 * $crlat2 + $this->b2 * $srlat2);
$this->ca = ($this->alt + $this->d) / $this->r;
$this->sa = $this->c2 * $crlat * $srlat / ($this->r * $this->d);
}
if ($this->glon != $this->olon){
for ($m = 2; $m <= $this->maxord; $m++){
$this->sp[$m] = $this->sp[1] * $this->cp[$m-1] + $this->cp[1] * $this->sp[$m-1];
$this->cp[$m] = $this->cp[1] * $this->cp[$m-1] - $this->sp[1] * $this->sp[$m-1];
}
}
$aor = $this->re / $this->r;
$ar = $aor * $aor;
$br = 0; $bt = 0; $bp = 0; $bpp = 0;
for($n = 1; $n <= $this->maxord; $n++){
$ar = $ar * $aor;
for ($m = 0, $D3 = 1, $D4 = ($n + $m + $D3) / $D3; $D4 > 0; $D4--, $m += $D3){
//COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
//AND DERIVATIVES VIA RECURSION RELATIONS
if($this->alt != $this->oalt || $this->glat != $this->olat){
if($n == $m){
$this->snorm[$n + $m * 13] = $this->st * $this->snorm[$n - 1 + ($m - 1) * 13];
$this->dp[$m][$n] = $this->st * $this->dp[$m-1][$n-1]+ $this->ct* $this->snorm[$n - 1 + ($m - 1) * 13];
}
if($n == 1 && $m == 0){
$this->snorm[$n + $m * 13] = $this->ct * $this->snorm[$n - 1 + $m * 13];
$this->dp[$m][$n] = $this->ct * $this->dp[$m][$n - 1] - $this->st * $this->snorm[$n - 1 + $m * 13];
}
if($n > 1 && $n != $m){
if($m > $n - 2)
$this->snorm[$n - 2 + $m * 13] = 0.0;
if($m > $n - 2)
$this->dp[$m][$n - 2] = 0.0;
$this->snorm[$n + $m * 13] = $this->ct * $this->snorm[$n - 1 + $m * 13] - $this->k[$m][$n] * $this->snorm[$n - 2 + $m * 13];
$this->dp[$m][$n] = $this->ct * $this->dp[$m][$n - 1] - $this->st * $this->snorm[$n - 1 + $m * 13] - $this->k[$m][$n] * $this->dp[$m][$n - 2];
}
}
//TIME ADJUST THE GAUSS COEFFICIENTS
if($this->time != $this->otime){
$this->tc[$m][$n] = $this->c[$m][$n] + $dt * $this->cd[$m][$n];
if($m != 0)
$this->tc[$n][$m - 1] = $this->c[$n][$m - 1]+ $dt * $this->cd[$n][$m - 1];
}
//ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
// $temp1=0; $temp2=0;
$par = $ar * $this->snorm[ $n + $m * 13];
if($m == 0){
$temp1 = $this->tc[$m][$n] * $this->cp[$m];
$temp2 = $this->tc[$m][$n] * $this->sp[$m];
} else{
$temp1 = $this->tc[$m][$n] * $this->cp[$m] + $this->tc[$n][$m - 1] * $this->sp[$m];
$temp2 = $this->tc[$m][$n] * $this->sp[$m] - $this->tc[$n][$m - 1] * $this->cp[$m];
}
$bt = $bt - $ar * $temp1 * $this->dp[$m][$n];
$bp += ($this->fm[$m] * $temp2 * $par);
$br += ($this->fn[$n] * $temp1 * $par);
//SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES
if($this->st == 0.0 && $m == 1){
if($n == 1)
$this->pp[n] = $this->pp[n - 1];
else
$this->pp[n] = $this->ct * $this->pp[n - 1] - $this->k[m][n] * $this->pp[n - 2];
$parp = $ar * $this->pp[$n];
$bpp += ($this->fm[$m] * $temp2 * $parp);
}
} //for(m...)
} //for(n...)
if($this->st == 0.0)
$bp = $bpp;
else
$bp /= $this->st;
//ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO
//GEODETIC COORDINATES
// bx must be the east-west field component
// by must be the north-south field component
// bz must be the vertical field component.
$this->bx = -$bt * $this->ca - $br * $this->sa;
$this->by = $bp;
$this->bz = $bt * $this->sa - $br * $this->ca;
//COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND
//TOTAL INTENSITY (TI)
$this->bh = sqrt(($this->bx * $this->bx)+($this->by * $this->by));
$this->ti = sqrt(($this->bh * $this->bh)+($this->bz * $this->bz));
// Calculate the declination.
$this->dec = (atan2($this->by, $this->bx) / $dtr);
printf ("Dec is: %f\n", $this->dec );
$this->dip = (atan2($this->bz, $this->bh) / $dtr);
// This is the variation for grid navigation.
// Not used at this time. See St. Ledger for explanation.
//COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT
//GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC
//(I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES)
// Grid North is referenced to the 0 Meridian of a polar
// stereographic projection.
//OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0
/*
gv = -999.0;
if (Math.abs(glat) >= 55.){
if (glat > 0.0 && glon >= 0.0)
gv = dec-glon;
if (glat > 0.0 && glon < 0.0)
gv = dec + Math.abs(glon);
if (glat < 0.0 && glon >= 0.0)
gv = dec+glon;
if (glat < 0.0 && glon < 0.0)
gv = dec - Math.abs(glon);
if (gv > +180.0)
gv -= 360.0;
if (gv < -180.0)
gv += 360.0;
}
*/
$this->otime = $this->time;
$this->oalt = $this->alt;
$this->olat = $this->glat;
$this->olon = $this->glon;
}
/**
* Returns the declination from the Department of
* Defense geomagnetic model and data, in degrees. The
* magnetic heading + declination = true heading.
* (True heading + variation = magnetic heading.)
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double All parameters together.
*/
public function getGeoMag($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return array(
'declination' => $this->dec,
'dip_angle' => $this->dip,
'intensity' => $this->ti,
'horizontal_intensity' => $this->bh,
'vertical_intensity' => $this->bz,
'north_intensity' => $this->bx,
'east_intensity' => $this->by,
);
}
/**
* Returns the declination from the Department of
* Defense geomagnetic model and data, in degrees. The
* magnetic heading + declination = true heading.
* (True heading + variation = magnetic heading.)
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double The declination in degrees.
*/
public function getDeclination($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
// printf("getDeclination(%f, %f) -> %f\n",$dlat,$dlong, $this->dec);
return $this->dec;
}
/**
* Returns the magnetic field intensity from the
* Department of Defense geomagnetic model and data
* in nano Tesla.
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double Magnetic field strength in nano Tesla.
*/
public function getIntensity($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return $this->ti;
}
/**
* Returns the horizontal magnetic field intensity from the
* Department of Defense geomagnetic model and data
* in nano Tesla.
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double The horizontal magnetic field strength in nano Tesla.
*/
public function getHorizontalIntensity($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return $this->bh;
}
/**
* Returns the vertical magnetic field intensity from the
* Department of Defense geomagnetic model and data
* in nano Tesla.
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double The vertical magnetic field strength in nano Tesla.
*/
public function getVerticalIntensity($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return $this->bz;
}
/**
* Returns the northerly magnetic field intensity from the
* Department of Defense geomagnetic model and data
* in nano Tesla.
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double The northerly component of the magnetic field strength in nano Tesla.
*/
public function getNorthIntensity($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return $this->bx;
}
/**
* Returns the easterly magnetic field intensity from the
* Department of Defense geomagnetic model and data
* in nano Tesla.
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double The easterly component of the magnetic field strength in nano Tesla.
*/
public function getEastIntensity($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return $this->by;
}
/**
* Returns the magnetic field dip angle from the
* Department of Defense geomagnetic model and data,
* in degrees.
*
* @param double dlat Latitude in decimal degrees.
* @param double dlong Longitude in decimal degrees.
* @param double year The date as a decimial year (optional)
* @param double altitude The altitude in kilometers (optional).
*
* @return double The magnetic field dip angle, in degrees.
*/
public function getDipAngle($dlat, $dlong, $year=null, $altitude=null)
{
$this->calcGeoMag($dlat, $dlong, $year, $altitude );
return $this->dip;
}
/** This method sets the input data to the internal fit coefficents.
* If there is an exception reading the input file WMM.COF, these values
* are used.
*
* NOTE: This method is not tested by the JUnit test, unless the WMM.COF file
* is missing.
*/
private function setCoeff()
{
$this->c[0][0] = 0.0;
$this->cd[0][0] = 0.0;
// $this->epoch = Double.parseDouble($this->input[0].trim().split("[\\s]+")[0]);
$this->epoch = doubleval(preg_split('/\s+/', " ".$this->input[0])[1]);
$this->defaultDate = $this->epoch + 2.5;
// String [] tokens;
//loop to get data from internal values
for($i=1; $i < count($this->input); $i++)
{
//tokens = input[i].trim().split("[\\s]+");
$tokens = preg_split('/\s+/', " ".$this->input[$i]);
$n = intval($tokens[1]);
$m = intval($tokens[2]);
$gnm = doubleval($tokens[3]);
$hnm = doubleval($tokens[4]);
$dgnm = doubleval($tokens[5]);
$dhnm = doubleval($tokens[6]);
if ($m <= $n)
{
$this->c [$m][$n] = $gnm;
$this->cd[$m][$n] = $dgnm;
if ($m != 0)
{
$this->c [$n][$m-1] = $hnm;
$this->cd[$n][$m-1] = $dhnm;
}
}
}
}
/**<p>
* Given a Gregorian Calendar object, this returns the decimal year
* value for the calendar, accurate to the day of the input calendar.
* The hours, minutes, and seconds of the date are ignored.</p><p>
*
* If the input Gregorian Calendar is new GregorianCalendar(2012, 6, 1), all of
* the first of July is counted, and this returns 2012.5. (183 days out of 366)</p><p>
*
* If the input Gregorian Calendar is new GregorianCalendar(2010, 0, 0), the first
* of January is not counted, and this returns 2010.0</p><p>
*
* @param int seconds (cal Has the date (year, month, and day of the month))
* @return int The date in decimal years
*/
// public double decimalYear(GregorianCalendar cal)
public function decimalYear($timestamp=null)
{
if ($timestamp === null) $timestamp = time();
$year = intval(date('Y',$timestamp)); // cal.get(Calendar.YEAR);
// double daysInYear;
if(intval(date('L',$timestamp))) {// (cal.isLeapYear(year))
$daysInYear = 366.0;
} else {
$daysInYear = 365.0;
}
return $year + intval(date('z',$timestamp))/ $daysInYear; //(cal.get(Calendar.DAY_OF_YEAR))/daysInYear;
}
}
?>
\ No newline at end of file
<?php
require_once ('geomag.php');
/**<p>
* Test values from
* <a href ="http://www.ngdc.noaa.gov/geomag/WMM/soft.shtml"> the National GeoPhysical Data Center.</a>.
* Click on the WMM2015testvalues.pdf link.
* </p><p>
* You have to run this test twice. Once with the WMM.COF
* file present, and then with it missing. Otherwise the
* setCoeff method is never tested.</p>
*
* @version 1.0 Apr 14, 2006
* @author John St. Ledger
* @version 1.1 Jan 28, 2009
* <p>Added 2006 test values.</p>
* @version 1.2 Jan 5, 2010
* <p>Updated with the test values for the 2010 WMM.COF coefficients. From page 18 of
* <i>The US/UK World Magnetic Model for 2010-2015, Preliminary online version containing
* final WMM2010 model coefficients</i></p>
* @version 1.3 Jan 15, 2015
* <p>Updated with the test values for the 2015 WMM.COF coefficients. From the test values WMM2015testvalues.pdf
* from the WMM web site.</p> * @version 1.4 May 26, 2015
* <p>Fixed the East-West, North-South bug discovered by Martin Frassl.</p>
*/
class TSAGeoMagTest
{
public function __construct(){
$this->failure = false;
$this->magModel = new GeoMag();
}
public function all_tests(){
$results = array();
$results['getDeclination'] = $this->testDeclination()?"pass":"fail";
$results['getDipAngle'] = $this->testDipAngle()?"pass":"fail";
$results['getHorizontalIntensity'] = $this->testHorizontalIntensity()?"pass":"fail";
$results['getNorthIntensity'] = $this->testNorthIntensity()?"pass":"fail";
$results['getEastIntensity'] = $this->testEastIntensity()?"pass":"fail";
$results['getVerticalIntensity'] = $this->testVerticalIntensity()?"pass":"fail";
$results['getIntensity'] = $this->testIntensity()?"pass":"fail";
$results['decimalYear'] = $this->testdecimalYear()?"pass":"fail";
print_r($results);
$all = $this->magModel->getGeoMag(80, 0, 2015, 0);
print_r($all);
$all = $this->magModel->getGeoMag(0, 120, 2015, 0);
print_r($all);
$all = $this->magModel->getGeoMag(40.723458, -111.932889, 2015, 0);
print_r($all);
$all = $this->magModel->getGeoMag(40.723458, -111.932889, 2020, 0);
print_r($all);
$all = $this->magModel->getGeoMag(40.723458, -111.932889, 2020, 1.5);
print_r($all);
return $results;
}
/**
* Test method for {@link d3.env.TSAGeoMag#getDeclination(double, double, double, double)}
*/
public function testDeclination()
{
$this->newTest();
$this->assertEquals(-3.85, $this->magModel->getDeclination(80, 0, 2015, 0) , 5.0E-03);
$this->assertEquals(0.57, $this->magModel->getDeclination(0, 120, 2015, 0) , 5.0E-03);
$this->assertEquals(69.81, $this->magModel->getDeclination(-80, 240, 2015, 0) , 5.0E-03);
$this->assertEquals(-4.27, $this->magModel->getDeclination(80, 0, 2015, 100) , 5.0E-03);
$this->assertEquals(0.56, $this->magModel->getDeclination(0, 120, 2015, 100) , 5.0E-03);
$this->assertEquals(69.22, $this->magModel->getDeclination(-80, 240, 2015, 100) , 5.0E-03);
$this->assertEquals(-2.75, $this->magModel->getDeclination(80, 0, 2017.5, 0) , 5.0E-03);
$this->assertEquals(0.32, $this->magModel->getDeclination(0, 120, 2017.5, 0) , 5.0E-03);
$this->assertEquals(69.58, $this->magModel->getDeclination(-80, 240, 2017.5, 0) , 5.0E-03);
$this->assertEquals(-3.17, $this->magModel->getDeclination(80, 0, 2017.5, 100) , 5.0E-03);
$this->assertEquals(0.32, $this->magModel->getDeclination(0, 120, 2017.5, 100) , 5.0E-03);
$this->assertEquals(69.00, $this->magModel->getDeclination(-80, 240, 2017.5, 100) , 5.0E-03);
$this->assertEquals(-2.75, $this->magModel->getDeclination(80, 0) , 5.0E-03);
$this->assertEquals(0.32, $this->magModel->getDeclination(0, 120) , 5.0E-03);
$this->assertEquals(69.58, $this->magModel->getDeclination(-80, 240) , 5.0E-03);
return $this->getSuccess();
}
/**
* Test method for {@link d3.env.TSAGeoMag#getDipAngle(double, double, double, double)}
*/
public function testDipAngle()
{
$this->newTest();
$this->assertEquals(83.04, $this->magModel->getDipAngle(80, 0, 2015, 0) , 5E-03);
$this->assertEquals(-15.89, $this->magModel->getDipAngle(0, 120, 2015, 0) , 5.0E-03);
$this->assertEquals(-72.39, $this->magModel->getDipAngle(-80, 240, 2015, 0) , 5.0E-03);
$this->assertEquals(83.09, $this->magModel->getDipAngle(80, 0, 2015, 100) , 5.0E-03);
$this->assertEquals(-16.01, $this->magModel->getDipAngle(0, 120, 2015, 100) , 5.0E-03);
$this->assertEquals(-72.57, $this->magModel->getDipAngle(-80, 240, 2015, 100) , 5.0E-03);
$this->assertEquals(83.08, $this->magModel->getDipAngle(80, 0, 2017.5, 0) , 5.0E-03);
$this->assertEquals(-15.57, $this->magModel->getDipAngle(0, 120, 2017.5, 0) , 5.0E-03);
$this->assertEquals(-72.28, $this->magModel->getDipAngle(-80, 240, 2017.5, 0) , 5.0E-03);
$this->assertEquals(83.13, $this->magModel->getDipAngle(80, 0, 2017.5, 100) , 5.0E-03);
$this->assertEquals(-15.70, $this->magModel->getDipAngle(0, 120, 2017.5, 100) , 5.0E-03);
$this->assertEquals(-72.45, $this->magModel->getDipAngle(-80, 240, 2017.5, 100) , 5.0E-03);
$this->assertEquals(83.08, $this->magModel->getDipAngle(80, 0) , 5.0E-03);
$this->assertEquals(-15.57, $this->magModel->getDipAngle(0, 120) , 5.0E-03);
$this->assertEquals(-72.28, $this->magModel->getDipAngle(-80, 240) , 5.0E-03);
return $this->getSuccess();
}
/**
* Test method for {@link d3.env.TSAGeoMag#getHorizontalIntensity(double, double, double, double)} in nT
* and {@link d3.env.TSAGeoMag#getHorizontalIntensity(double, double)} in nT
*/
public function testHorizontalIntensity()
{
$this->newTest();
$this->assertEquals(6642.1, $this->magModel->getHorizontalIntensity(80, 0, 2015, 0) , 5.0E-02);
$this->assertEquals(39520.2, $this->magModel->getHorizontalIntensity(0, 120, 2015, 0) , 5.0E-02);
$this->assertEquals(16793.5, $this->magModel->getHorizontalIntensity(-80, 240, 2015, 0) , 5.0E-02);
$this->assertEquals(6331.9, $this->magModel->getHorizontalIntensity(80, 0, 2015, 100) , 5.0E-02);
$this->assertEquals(37537.3, $this->magModel->getHorizontalIntensity(0, 120, 2015, 100) , 5.0E-02);
$this->assertEquals(15820.7, $this->magModel->getHorizontalIntensity(-80, 240, 2015, 100) , 5.0E-02);
$this->assertEquals(6607.0, $this->magModel->getHorizontalIntensity(80, 0, 2017.5, 0) , 5.0E-02);
$this->assertEquals(39572.0, $this->magModel->getHorizontalIntensity(0, 120, 2017.5, 0) , 5.0E-02);
$this->assertEquals(16839.1, $this->magModel->getHorizontalIntensity(-80, 240, 2017.5, 0) , 5.0E-02);
$this->assertEquals(6300.1, $this->magModel->getHorizontalIntensity(80, 0, 2017.5, 100) , 5.0E-02);
$this->assertEquals(37586.1, $this->magModel->getHorizontalIntensity(0, 120, 2017.5, 100) , 5.0E-02);
$this->assertEquals(15862.0, $this->magModel->getHorizontalIntensity(-80, 240, 2017.5, 100) , 5.0E-02);
$this->assertEquals(6607.0, $this->magModel->getHorizontalIntensity(80, 0) , 5.0E-02);
$this->assertEquals(39572.0, $this->magModel->getHorizontalIntensity(0, 120) , 5.0E-02);
$this->assertEquals(16839.1, $this->magModel->getHorizontalIntensity(-80, 240) , 5.0E-02);
return $this->getSuccess();
}
/**
* Test method for d3.env.TSAGeoMag.getEastIntensity() in nT
*/
public function testNorthIntensity()
{
$this->newTest();
$this->assertEquals(6627.1, $this->magModel->getNorthIntensity(80, 0, 2015, 0) , 5.0E-02);
$this->assertEquals(39518.2, $this->magModel->getNorthIntensity(0, 120, 2015, 0) , 5.0E-02);
$this->assertEquals(5797.3, $this->magModel->getNorthIntensity(-80, 240, 2015, 0) , 5.0E-02);
$this->assertEquals(6314.3, $this->magModel->getNorthIntensity(80, 0, 2015, 100) , 5.0E-02);
$this->assertEquals(37535.6, $this->magModel->getNorthIntensity(0, 120, 2015, 100) , 5.0E-02);
$this->assertEquals(5613.1, $this->magModel->getNorthIntensity(-80, 240, 2015, 100) , 5.0E-02);
$this->assertEquals(6599.4, $this->magModel->getNorthIntensity(80, 0, 2017.5, 0) , 5.0E-02);
$this->assertEquals(39571.4, $this->magModel->getNorthIntensity(0, 120, 2017.5, 0) , 5.0E-02);
$this->assertEquals(5873.8, $this->magModel->getNorthIntensity(-80, 240, 2017.5, 0) , 5.0E-02);
$this->assertEquals(6290.5, $this->magModel->getNorthIntensity(80, 0, 2017.5, 100) , 5.0E-02);
$this->assertEquals(37585.5, $this->magModel->getNorthIntensity(0, 120, 2017.5, 100) , 5.0E-02);
$this->assertEquals(5683.5, $this->magModel->getNorthIntensity(-80, 240, 2017.5, 100) , 5.0E-02);
$this->assertEquals(6599.4, $this->magModel->getNorthIntensity(80, 0) , 5.0E-02);
$this->assertEquals(39571.4, $this->magModel->getNorthIntensity(0, 120) , 5.0E-02);
$this->assertEquals(5873.8, $this->magModel->getNorthIntensity(-80, 240) , 5.0E-02);
return $this->getSuccess();
}
/**
* Test method for d3.env.TSAGeoMag.getNorthIntensity() in nT
*/
public function testEastIntensity()
{
$this->newTest();
$this->assertEquals( -445.9, $this->magModel->getEastIntensity(80, 0, 2015, 0) , 5.0E-02);
$this->assertEquals( 392.9, $this->magModel->getEastIntensity(0, 120, 2015, 0) , 5.0E-02);
$this->assertEquals(15761.1, $this->magModel->getEastIntensity(-80, 240, 2015, 0) , 5.0E-02);
$this->assertEquals( -471.6, $this->magModel->getEastIntensity(80, 0, 2015, 100) , 5.0E-02);
$this->assertEquals( 364.4, $this->magModel->getEastIntensity(0, 120, 2015, 100) , 5.0E-02);
$this->assertEquals(14791.5, $this->magModel->getEastIntensity(-80, 240, 2015, 100) , 5.0E-02);
$this->assertEquals( -317.1, $this->magModel->getEastIntensity(80, 0, 2017.5, 0) , 5.0E-02);
$this->assertEquals( 222.5, $this->magModel->getEastIntensity(0, 120, 2017.5, 0) , 5.0E-02);
$this->assertEquals(15781.4, $this->magModel->getEastIntensity(-80, 240, 2017.5, 0) , 5.0E-02);
$this->assertEquals( -348.5, $this->magModel->getEastIntensity(80, 0, 2017.5, 100) , 5.0E-02);
$this->assertEquals( 209.5, $this->magModel->getEastIntensity(0, 120, 2017.5, 100) , 5.0E-02);
$this->assertEquals(14808.8, $this->magModel->getEastIntensity(-80, 240, 2017.5, 100) , 5.0E-02);
$this->assertEquals( -317.1, $this->magModel->getEastIntensity(80, 0) , 5.0E-02);
$this->assertEquals( 222.5, $this->magModel->getEastIntensity(0, 120) , 5.0E-02);
$this->assertEquals(15781.4, $this->magModel->getEastIntensity(-80, 240) , 5.0E-02);
return $this->getSuccess();
}
/**
* Test method for d3.env.TSAGeoMag.getVerticalIntensity()
*/
public function testVerticalIntensity()
{
$this->newTest();
$this->assertEquals( 54432.3, $this->magModel->getVerticalIntensity(80, 0, 2015, 0) , 5.0E-02);
$this->assertEquals(-11252.4, $this->magModel->getVerticalIntensity(0, 120, 2015, 0) , 5.0E-02);
$this->assertEquals(-52919.1, $this->magModel->getVerticalIntensity(-80, 240, 2015, 0) , 5.0E-02);
$this->assertEquals( 52269.8, $this->magModel->getVerticalIntensity(80, 0, 2015, 100) , 5.0E-02);
$this->assertEquals(-10773.4, $this->magModel->getVerticalIntensity(0, 120, 2015, 100) , 5.0E-02);
$this->assertEquals(-50378.6, $this->magModel->getVerticalIntensity(-80, 240, 2015, 100) , 5.0E-02);
$this->assertEquals( 54459.2, $this->magModel->getVerticalIntensity(80, 0, 2017.5, 0) , 5.0E-02);
$this->assertEquals(-11030.1, $this->magModel->getVerticalIntensity(0, 120, 2017.5, 0) , 5.0E-02);
$this->assertEquals(-52687.9, $this->magModel->getVerticalIntensity(-80, 240, 2017.5, 0) , 5.0E-02);
$this->assertEquals( 52292.7, $this->magModel->getVerticalIntensity(80, 0, 2017.5, 100) , 5.0E-02);
$this->assertEquals(-10564.2, $this->magModel->getVerticalIntensity(0, 120, 2017.5, 100) , 5.0E-02);
$this->assertEquals(-50163.0, $this->magModel->getVerticalIntensity(-80, 240, 2017.5, 100) , 5.0E-02);
$this->assertEquals( 54459.2, $this->magModel->getVerticalIntensity(80, 0) , 5.0E-02);
$this->assertEquals(-11030.1, $this->magModel->getVerticalIntensity(0, 120) , 5.0E-02);
$this->assertEquals(-52687.9, $this->magModel->getVerticalIntensity(-80, 240) , 5.0E-02);
return $this->getSuccess();
}
/**
* Test method for d3.env.TSAGeoMag.getIntensity()
*/
public function testIntensity()
{
$this->newTest();
$this->assertEquals(54836.0, $this->magModel->getIntensity(80, 0, 2015, 0) , 5.0E-02);
$this->assertEquals(41090.9, $this->magModel->getIntensity(0, 120, 2015, 0) , 5.0E-02);
$this->assertEquals(55519.8, $this->magModel->getIntensity(-80, 240, 2015, 0) , 5.0E-02);
$this->assertEquals(52652.0, $this->magModel->getIntensity(80, 0, 2015, 100) , 5.0E-02);
$this->assertEquals(39052.7, $this->magModel->getIntensity(0, 120, 2015, 100) , 5.0E-02);
$this->assertEquals(52804.4, $this->magModel->getIntensity(-80, 240, 2015, 100) , 5.0E-02);
$this->assertEquals(54858.5, $this->magModel->getIntensity(80, 0, 2017.5, 0) , 5.0E-02);
$this->assertEquals(41080.5, $this->magModel->getIntensity(0, 120, 2017.5, 0) , 5.0E-02);
$this->assertEquals(55313.4, $this->magModel->getIntensity(-80, 240, 2017.5, 0) , 5.0E-02);
$this->assertEquals(52670.9, $this->magModel->getIntensity(80, 0, 2017.5, 100) , 5.0E-02);
$this->assertEquals(39042.5, $this->magModel->getIntensity(0, 120, 2017.5, 100) , 5.0E-02);
$this->assertEquals(52611.1, $this->magModel->getIntensity(-80, 240, 2017.5, 100) , 5.0E-02);
$this->assertEquals(54858.5, $this->magModel->getIntensity(80, 0) , 5.0E-02);
$this->assertEquals(41080.5, $this->magModel->getIntensity(0, 120) , 5.0E-02);
$this->assertEquals(55313.4, $this->magModel->getIntensity(-80, 240) , 5.0E-02);
//assertEquals(52672.8, $this->magModel->getIntensity(40, -105, 2014, 0), 0.05);
//assertEquals(52672.8, $this->magModel->getIntensity(40, -105, $this->magModel->decimalYear(new GregorianCalendar(2014, 0, 0)), 0), 0.05);
//assertEquals(52672.5, $this->magModel->getIntensity(40, -105, $this->magModel->decimalYear(new GregorianCalendar(2014, 0, 1)), 0), 0.05);
return $this->getSuccess();
}
/**
* test method for {@link d3.env.TSAGeoMag#decimalYear(GregorianCalendar)}
*/
public function testdecimalYear()
{
$this->newTest();
$mag = new GeoMag();
// GregorianCalendar cal = new GregorianCalendar(2010, 0, 0);
// assertEquals(2010.0, $mag->decimalYear(cal), 0.0);
$ts = strtotime("2010/1/1");
$this->assertEquals(2010.0, $mag->decimalYear($ts), 0.0);
// GregorianCalendar cal2 = new GregorianCalendar(2012, 6, 1); // the full day of July 1, 0 hours into 2 July
// assertTrue(cal2.isLeapYear(2012));
// assertEquals(2012.5, $mag->decimalYear(cal2), 0.0);
$ts = strtotime("2012/7/2");
$this->assertEquals(2012.5, $mag->decimalYear($ts), 0.0005);
// cal2 = new GregorianCalendar(2013, 3, 13);
// assertFalse(cal2.isLeapYear(2013));
// assertEquals(2013.282, $mag->decimalYear(cal2), 0.0005);
$ts = strtotime("2013/4/14");
$this->assertEquals(2013.282, $mag->decimalYear($ts), 0.0005);
return $this->getSuccess();
}
private function assertEquals ($a, $b, $tolerance) {
$ok = abs($a - $b) <= $tolerance;
$this->failure |= !$ok;
printf("assertEquals (%f, %f, %f) -> %d (%d)\n",$a, $b, $tolerance, $ok, $this->failure);
return $ok;
}
private function assertFalse ($f) {
$this->failure &= $f;
return !$f;
}
private function assertTrue ($t) {
$this->failure &= !$t;
return $t;
}
private function newTest(){
$this->failure = false;
}
private function getSuccess(){
return !$this->failure;
}
}
$tests = new TSAGeoMagTest();
$tests->all_tests();
?>
<?php
$phone_ip= "http://192.168.7.222";
$phone_port= "8123";
$phone_url = $phone_ip.":".$phone_port;
$log_file = fopen("/home/eyesis/git/motosat/attic/logs/phone07.log","w");
$ts_l = 14;
$acc_l = 35;
$mag_l = 35;
$long_l = 12;
$lat_l = 12;
$salign = 128; // used 108
//$ll = $ts_l + $acc_l + $mag_l + $long_l + $lat_l;
$fs = "%".$ts_l."s,%".$acc_l."s,%".$mag_l."s,%".$lat_l."s,%".$long_l."s";
libxml_use_internal_errors(true); // supress XML errors
while (true){
$xml = simplexml_load_file($phone_url);
//PHP Warning: simplexml_load_file(http://192.168.7.222:8123): failed to open stream: Connection refused in /home/eyesis/git/motosat/log_phone_sensors.php on line 18
//PHP Warning: simplexml_load_file(): I/O warning : failed to load external entity "http://192.168.7.222:8123" in /home/eyesis/git/motosat/log_phone_sensors.php on line 18
if ($xml) {
$ts = $xml->timestmap;
$acc = $xml->accelerometer;
$mag = $xml->magnetometer;
$long = $xml->longitude;
$lat = $xml->latitude;
$line = sprintf($fs, $ts, $acc, $mag, $long, $lat);
$line .= sprintf("%" . ($salign - strlen($line) - 1) . "s", " ") . "\n";
fwrite($log_file, $line);
}
}
function getDishAngles($xml_state){
$DishAngles = explode(",",$xml_state->DishAngle);
foreach ($DishAngles as $key => $value){
$DishAngles[$key] = floatval($value);
}
return $DishAngles;
}
/*
$xml=simplexml_load_string("<sensors><timestmap>1570632414615</timestmap><accelerometer>0.095215, -0.811707, 9.966614</accelerometer><magnetometer>-15.299988, 4.798889, 46.078491</magnetometer><azimuth>61.597424</azimuth><pitch>4.655828</pitch><roll>-0.547352</roll><longitude>-111.932889</longitude><latitude>40.723488</latitude><date>Wed Oct 09 14:46:54 MDT 2019</date></sensors>");
*/
//1570632414615 // 13
//-111.932889 // 11
//0.095215, -0.811707, 9.966614 // 29
// -15.299988, 4.798889, 46.078491 // 31
?>
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment