:return
         D is declination (+ve east)
         I is inclination (+ve down)
         H is horizontal intensity
         X is north component
         Y is east component
         Z is vertical component (+ve down)
         F is total intensity
556004
def igrf_value(lat, lon, alt=0., year=2005.):
    
    X, Y, Z, F = calculate.igrf12syn(0, year, 1, alt, lat, lon)
    D = FACT * np.arctan2(Y, X)
    H = np.sqrt(X * X + Y * Y)
    I = FACT * np.arctan2(Z, H)
    return D, I, H, X, Y, Z, F
509079
def angle2d(self):
        """"""
        if self.x==0:
            if   self.y<0:  return       math.pi/2.0*3
            elif self.y>0:  return       math.pi/2.0
            else:           return 0
        elif self.y==0:
            if   self.x<0:  return       math.pi
           #elif self.x>0:  return 0
            else:           return 0
        ans = math.atan( self.y / self.x )
        if self.x > 0:
            if self.y>0:    return ans
            else:           return ans + math.pi*2.0
        else:               return ans + math.pi
912894
def obfn_dfd(self):
        
        """

        Ef = sl.inner(self.Df, self.obfn_fvarf(), axis=self.cri.axisM) - \
            self.Sf
        return sl.rfl2norm2(Ef, self.S.shape, axis=self.cri.axisN) / 2.0
211962
function getDirectionOfLocation([x1, y1], [x2, y2]) {
  if (Math.abs(x2 - x1) > Math.abs(y2 - y1)) {
    if (x1 > x2) {
      return EAST;
    }

    return WEST;
  }

  if (y1 > y2) {
    return SOUTH;
  }

  return NORTH;
}
35836
def eqCoords(lon, lat):
    
    # Convert to radians
    _lambda = math.radians(lon)
    _beta = math.radians(lat)
    _epson = math.radians(23.44)  # The earth's inclination
    
    # Declination in radians
    decl = math.asin(math.sin(_epson) * math.sin(_lambda) * math.cos(_beta) + \
           math.cos(_epson) * math.sin(_beta))
    
    # Equatorial Distance in radians
    ED = math.acos(math.cos(_lambda) * math.cos(_beta) / math.cos(decl))
    
    # RA in radians
    ra = ED if lon < 180 else math.radians(360) - ED
    
    # Correctness of RA if longitude is close to 0º or 180º in a radius of 5º
    if (abs(angle.closestdistance(lon, 0)) < 5 or 
        abs(angle.closestdistance(lon, 180)) < 5):
            a = math.sin(ra) * math.cos(decl)
            b = math.cos(_epson) * math.sin(_lambda) * math.cos(_beta) - \
                math.sin(_epson) * math.sin(_beta)
            if (math.fabs(a-b) > 0.0003):
                ra = math.radians(360) - ra
                
    return (math.degrees(ra), math.degrees(decl))
859592
def dot(self, rhs):
    

    return self.x * rhs.x + self.y * rhs.y + self.z * rhs.z
1089007
def _derY(self,x,y):
        '''
        
        '''
        x_pos, y_pos = self.findSector(x,y)
        alpha, beta = self.findCoords(x,y,x_pos,y_pos)

        # Get four corners data for each point
        xA = self.x_values[x_pos,y_pos]
        xB = self.x_values[x_pos+1,y_pos]
        xC = self.x_values[x_pos,y_pos+1]
        xD = self.x_values[x_pos+1,y_pos+1]
        yA = self.y_values[x_pos,y_pos]
        yB = self.y_values[x_pos+1,y_pos]
        yC = self.y_values[x_pos,y_pos+1]
        yD = self.y_values[x_pos+1,y_pos+1]
        fA = self.f_values[x_pos,y_pos]
        fB = self.f_values[x_pos+1,y_pos]
        fC = self.f_values[x_pos,y_pos+1]
        fD = self.f_values[x_pos+1,y_pos+1]

        # Calculate components of the alpha,beta --> x,y delta translation matrix
        alpha_x = (1-beta)*(xB-xA) + beta*(xD-xC)
        alpha_y = (1-beta)*(yB-yA) + beta*(yD-yC)
        beta_x  = (1-alpha)*(xC-xA) + alpha*(xD-xB)
        beta_y  = (1-alpha)*(yC-yA) + alpha*(yD-yB)

        # Invert the delta translation matrix into x,y --> alpha,beta
        det = alpha_x*beta_y - beta_x*alpha_y
        y_alpha = -beta_x/det
        y_beta  = alpha_x/det

        # Calculate the derivative of f w.r.t. alpha and beta
        dfda = (1-beta)*(fB-fA) + beta*(fD-fC)
        dfdb = (1-alpha)*(fC-fA) + alpha*(fD-fB)

        # Calculate the derivative with respect to x (and return it)
        dfdy = y_alpha*dfda + y_beta*dfdb
        return dfdy
922125
function toLatLng()
	{
		$wgs84 = new ReferenceEllipsoid(ReferenceEllipsoid::WGS84_MAJ, ReferenceEllipsoid::WGS84_MIN);

		$UTM_F0	 = 0.9996;
		$a = $wgs84->maj;
		$eSquared = $wgs84->ecc;
		$ePrimeSquared = $eSquared / (1.0 - $eSquared);
		$e1 = (1 - sqrt(1 - $eSquared)) / (1 + sqrt(1 - $eSquared));

		$x = $this->easting - 500000.0;;
		$y = $this->northing;

		$zoneNumber = $this->lngZone;
		$zoneLetter = $this->latZone;

		$longitudeOrigin = ($zoneNumber - 1.0) * 6.0 - 180.0 + 3.0;

		// Correct y for southern hemisphere
		if ((ord($zoneLetter) - ord("N")) < 0) {
			$y -= 10000000.0;
		}

		$m = $y / $UTM_F0;

		$mu = $m
			/ ($a
			* (1.0
				- $eSquared / 4.0
				- 3.0 * $eSquared * $eSquared / 64.0
				- 5.0
				* pow($eSquared, 3.0)
				/ 256.0));

		$phi1Rad = $mu
			+ (3.0 * $e1 / 2.0 - 27.0 * pow($e1, 3.0) / 32.0) * sin(2.0 * $mu)
			+ (21.0 * $e1 * $e1 / 16.0 - 55.0 * pow($e1, 4.0) / 32.0)
			* sin(4.0 * $mu)
			+ (151.0 * pow($e1, 3.0) / 96.0) * sin(6.0 * $mu);

		$n = $a / sqrt(1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad));

		$t = tan($phi1Rad) * tan($phi1Rad);
		$c = $ePrimeSquared * cos($phi1Rad) * cos($phi1Rad);

		$r = $a
			* (1.0 - $eSquared)
			/ pow(
			1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad),
			1.5);

		$d = $x / ($n * $UTM_F0);

		$latitude = ($phi1Rad
			- ($n * tan($phi1Rad) / $r)
			* ($d * $d / 2.0
				- (5.0
				+ (3.0 * $t)
				+ (10.0 * $c)
				- (4.0 * $c * $c)
				- (9.0 * $ePrimeSquared))
				* pow($d, 4.0)
				/ 24.0
				+ (61.0
				+ (90.0 * $t)
				+ (298.0 * $c)
				+ (45.0 * $t * $t)
				- (252.0 * $ePrimeSquared)
				- (3.0 * $c * $c))
				* pow($d, 6.0)
				/ 720.0)) * (180.0 / pi());

		$longitude = $longitudeOrigin + (($d
			- (1.0 + 2.0 * $t + $c) * pow($d, 3.0) / 6.0
			+ (5.0
			- (2.0 * $c)
			+ (28.0 * $t)
			- (3.0 * $c * $c)
			+ (8.0 * $ePrimeSquared)
			+ (24.0 * $t * $t))
			* pow($d, 5.0)
			/ 120.0)
			/ cos($phi1Rad)) * (180.0 / pi());

		return new LatLng($latitude, $longitude);
	}
38496
def obfn_dfd(self):
        
        """

        DXBf = sl.dot(self.B, sl.inner(self.Df, self.obfn_fvarf(),
                                       axis=self.cri.axisM),
                       axis=self.cri.axisC)
        Ef = DXBf - self.Sf
        return sl.rfl2norm2(Ef, self.S.shape, axis=self.cri.axisN) / 2.0
855504
def dec2hms(dec):
    
    DEGREE = 360.
    HOUR = 24.
    MINUTE = 60.
    SECOND = 3600.
    
    dec = float(dec)
    fhour = dec*(HOUR/DEGREE)
    hour = int(fhour)

    fminute = (fhour - hour)*MINUTE
    minute = int(fminute)
    
    second = (fminute - minute)*MINUTE
    return (hour, minute, second)
190326
function(a, b, p, dir) {

      var s = this['SCALE'];
      a = mod(a, s);
      b = mod(b, s);

      if (a === b)
        return a;

      // dir becomes an offset if we have to add a full revolution (=scale)
      if (!dir)
        dir =-s;
      else if ((dir === 1) === (a < b))
        dir*= s;
      else
        dir = 0;

      return mod(a + p * (b - a - dir), s);
    }
348641
def _getDirection(coord1, coord2):
    
    x1, y1 = coord1
    x2, y2 = coord2

    if x1 == x2 and y1 == y2:
        return None # two coordinates are the same.
    elif x1 == x2 and y1 > y2:
        return UP
    elif x1 == x2 and y1 < y2:
        return DOWN
    elif x1 > x2 and y1 == y2:
        return LEFT
    elif x1 < x2 and y1 == y2:
        return RIGHT

    slope = float(y2 - y1) / float(x2 - x1)

    # Figure out which quadrant the line is going in, and then
    # determine the closest direction by calculating the slope
    if x2 > x1 and y2 < y1: # up right quadrant
        if slope > -0.4142:
            return RIGHT # slope is between 0 and 22.5 degrees
        elif slope < -2.4142:
            return UP # slope is between 67.5 and 90 degrees
        else:
            return UPRIGHT # slope is between 22.5 and 67.5 degrees
    elif x2 > x1 and y2 > y1: # down right quadrant
        if slope > 2.4142:
            return DOWN
        elif slope < 0.4142:
            return RIGHT
        else:
            return DOWNRIGHT
    elif x2 < x1 and y2 < y1: # up left quadrant
        if slope < 0.4142:
            return LEFT
        elif slope > 2.4142:
            return UP
        else:
            return UPLEFT
    elif x2 < x1 and y2 > y1: # down left quadrant
        if slope < -2.4142:
            return DOWN
        elif slope > -0.4142:
            return LEFT
        else:
            return DOWNLEFT
913110
def rotateInDeclination(v1, theta_deg):
    """"""
    axis = np.array([0,-1,0])
    return rotateAroundVector(v1, axis, theta_deg)
792102
def inclination(vector)
      Math.atan2(-vector.y, Math.sqrt(vector.x * vector.x + vector.z * vector.z))
    end
922246
function calcDirection(start, end) {
  var direction = new Vec2(end);
  direction.subtract(start);
  direction.normalize();
  return direction;
}
1041941
def dxdy(line):
    
    x0, y0, x1, y1 = line
    dx = float(x1 - x0)
    dy = float(y1 - y0)
    f = hypot(dx, dy)
    return dx / f, dy / f
1160836
def focus0(self):
        '''
        

        '''
        f = Point(self.center)

        if self.xAxisIsMajor:
            f.x -= self.linearEccentricity
        else:
            f.y -= self.linearEccentricity
        return f