Interpolate using cubic spline approximation
        Returns the same as interpolate(lons,lats,data,order=3)
505101
def interpolate_cubic(self, lons, lats, data):
        
        return self.interpolate(lons, lats, data, order=3)
1203603
private int getCubicInterpolated(final int currentSamplePos, final int currentTuningPos)
	{
		final int poslo = (currentTuningPos >> CubicSpline.SPLINE_FRACSHIFT) & CubicSpline.SPLINE_FRACMASK;
		
		final long v1 = (((currentSamplePos-1)<0)?0L:	(long)CubicSpline.lut[poslo  ]*(long)sample[currentSamplePos-1]) +
				  		(								(long)CubicSpline.lut[poslo+1]*(long)sample[currentSamplePos  ]) +
				  		(								(long)CubicSpline.lut[poslo+2]*(long)sample[currentSamplePos+1]) +
				  		(								(long)CubicSpline.lut[poslo+3]*(long)sample[currentSamplePos+2]);
		
		return (int)(v1 >> CubicSpline.SPLINE_QUANTBITS);
	}
401382
private double cubicInterpolate(final double Y0, final double Y1, final double Y2, final double Y3, final double MU) {
        final double A0;
        final double A1;
        final double A2;
        final double A3;
        final double MU2;

        MU2 = MU * MU;
        A0 = Y3 - Y2 - Y0 + Y1;
        A1 = Y0 - Y1 - A0;
        A2 = Y2 - Y0;
        A3 = Y1;

        return ((A0 * MU * MU2) + (A1 * MU2) + (A2 * MU) + A3);
    }
447951
function getEquallySpacedSlot(x, y, from, to, numberOfPoints) {
                var xLength = x.length;

                var step = (to - from) / (numberOfPoints - 1);
                var halfStep = step / 2;
                var lastStep = x[x.length - 1] - x[x.length - 2];

                var start = from - halfStep;
                var output = new Array(numberOfPoints);

                // Init main variables
                var min = start;
                var max = start + step;

                var previousX = -Number.MAX_VALUE;
                var previousY = 0;
                var nextX = x[0];
                var nextY = y[0];
                var frontOutsideSpectra = 0;
                var backOutsideSpectra = true;

                var currentValue = 0;

                // for slot algorithm
                var currentPoints = 0;

                var i = 1; // index of input
                var j = 0; // index of output

                main: while (true) {
                    if (previousX >= nextX) throw (new Error('x must be an increasing serie'));
                    while (previousX - max > 0) {
                        // no overlap with original point, just consume current value
                        if (backOutsideSpectra) {
                            currentPoints++;
                            backOutsideSpectra = false;
                        }

                        output[j] = currentPoints <= 0 ? 0 : currentValue / currentPoints;
                        j++;

                        if (j === numberOfPoints) {
                            break main;
                        }

                        min = max;
                        max += step;
                        currentValue = 0;
                        currentPoints = 0;
                    }

                    if (previousX > min) {
                        currentValue += previousY;
                        currentPoints++;
                    }

                    if (previousX === -Number.MAX_VALUE || frontOutsideSpectra > 1) {
                        currentPoints--;
                    }

                    previousX = nextX;
                    previousY = nextY;

                    if (i < xLength) {
                        nextX = x[i];
                        nextY = y[i];
                        i++;
                    } else {
                        nextX += lastStep;
                        nextY = 0;
                        frontOutsideSpectra++;
                    }
                }

                return output;
            }
1294724
def interpolate_xarray(xpoints, ypoints, values, shape, kind='cubic',
                       blocksize=CHUNK_SIZE):
    """"""
    vchunks = range(0, shape[0], blocksize)
    hchunks = range(0, shape[1], blocksize)

    token = tokenize(blocksize, xpoints, ypoints, values, kind, shape)
    name = 'interpolate-' + token

    from scipy.interpolate import interp2d
    interpolator = interp2d(xpoints, ypoints, values, kind=kind)

    dskx = {(name, i, j): (interpolate_slice,
                           slice(vcs, min(vcs + blocksize, shape[0])),
                           slice(hcs, min(hcs + blocksize, shape[1])),
                           interpolator)
            for i, vcs in enumerate(vchunks)
            for j, hcs in enumerate(hchunks)
            }

    res = da.Array(dskx, name, shape=list(shape),
                   chunks=(blocksize, blocksize),
                   dtype=values.dtype)
    return DataArray(res, dims=('y', 'x'))
1172933
function CubicBezier(coordinates) {
        if (typeof coordinates === "string") {
            this.coordinates = coordinates.split(",");
        } else {
            this.coordinates = coordinates;
        }

        if (!this.coordinates) {
            throw "No offsets were defined";
        }

        this.coordinates = this.coordinates.map(function (n) { return +n; });

        var i;
        for (i = 3; i >= 0; i--) {
            var xy = this.coordinates[i];
            if (isNaN(xy) || (((i % 2) === 0) && (xy < 0 || xy > 1))) {
                throw "Wrong coordinate at " + i + "(" + xy + ")";
            }
        }
    }
1319657
def lonlat2xyz(lons, lats):
    """"""
    R = 6370997.0
    x_coords = R * da.cos(da.deg2rad(lats)) * da.cos(da.deg2rad(lons))
    y_coords = R * da.cos(da.deg2rad(lats)) * da.sin(da.deg2rad(lons))
    z_coords = R * da.sin(da.deg2rad(lats))
    return x_coords, y_coords, z_coords
1097608
def to_lat_lon
      k0 = 0.9996
      x  = @e - 500000 # Remove Longitude offset
      y  = @n

      # Set hemisphere (1=Northern, 0=Southern)
      y    -= 10000000.0 unless UTMZones::northern_hemisphere? @zone

      ecc = @ellipsoid.eccentricity
      eccPrimeSquared = (ecc)/(1-ecc)
      m  = y / k0
      mu = m / (@ellipsoid.radius * (1 - ecc / 4 - 3 * ecc ** 2 / 64 - 5 * ecc ** 3 / 256))
      e1 = (1 - Math::sqrt(1 - ecc)) / (1 + Math::sqrt(1 - ecc))
      phi1rad = mu +
        (3 * e1 / 2 - 27 * e1 ** 3 / 32) * Math::sin(2 * mu) +
        (21 * e1 ** 2 / 16 - 55 * e1 ** 4 / 32) * Math::sin(4 * mu) +
        (151 * e1 ** 3 / 96) * Math::sin(6 * mu)
      n1 = @ellipsoid.radius / Math::sqrt(1 - ecc * Math::sin(phi1rad) ** 2)
      t1 = Math::tan(phi1rad) ** 2
      c1 = ecc * Math::cos(phi1rad)**2
      r1 = @ellipsoid.radius * (1 - ecc) / ((1 - ecc * Math::sin(phi1rad) ** 2) ** 1.5)
      d = x / (n1 * k0)
      latitude_rad = phi1rad-
        (n1 * Math::tan(phi1rad) / r1) * (d * d / 2 - 
                                    (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * eccPrimeSquared) * d ** 4 / 24 +
                                    (61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 
                                     252 * eccPrimeSquared - 3 * c1 * c1) * d ** 6 / 720)
      latitude_deg = latitude_rad * 180.0 / Math::PI
      lon_tmp = (d - (1 + 2 * t1 + c1) * d ** 3 / 6 +
                   (5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * eccPrimeSquared +
                    24 * t1 * t1) * d ** 5 / 120) / Math::cos(phi1rad)
      longitude_deg = UTMZones::lon_origin(@zone) + lon_tmp * 180.0 / Math::PI
      LatLon.new latitude_deg, longitude_deg
    end
638433
function getRandomPointOnSphere(r) {
  var u = THREE.Math.randFloat(0, 1);
  var v = THREE.Math.randFloat(0, 1);
  var theta = 2 * Math.PI * u;
  var phi = Math.acos(2 * v - 1);
  var x = r * Math.sin(theta) * Math.sin(phi);
  var y = r * Math.cos(theta) * Math.sin(phi);
  var z = r * Math.cos(phi);
  return {
    x,
    y,
    z,
  };
}
564172
def eqCoords(self, zerolat=False):
        
        lat = 0.0 if zerolat else self.lat
        return utils.eqCoords(self.lon, lat)
607014
def getAltitudeFromLatLon(self, lat, lon):
        
        # print "-----\nFromLatLon", lon, lat
        lat -= self.lat
        lon -= self.lon
        # print "lon, lat", lon, lat
        if lat < 0.0 or lat >= 1.0 or lon < 0.0 or lon >= 1.0:
            raise WrongTileError(self.lat, self.lon, self.lat+lat, self.lon+lon)
        x = lon * (self.size - 1)
        y = lat * (self.size - 1)
        # print "x,y", x, y
        x_int = int(x)
        x_frac = x - int(x)
        y_int = int(y)
        y_frac = y - int(y)
        # print "frac", x_int, x_frac, y_int, y_frac
        value00 = self.getPixelValue(x_int, y_int)
        value10 = self.getPixelValue(x_int+1, y_int)
        value01 = self.getPixelValue(x_int, y_int+1)
        value11 = self.getPixelValue(x_int+1, y_int+1)
        value1 = self._avg(value00, value10, x_frac)
        value2 = self._avg(value01, value11, x_frac)
        value  = self._avg(value1,  value2, y_frac)
        # print "%4d %4d | %4d\n%4d %4d | %4d\n-------------\n%4d" % (
        #        value00, value10, value1, value01, value11, value2, value)
        return value
551699
def interpolate(self, doy, depth, lat, lon, var):
        
        subset, dims = self.crop(doy, depth, lat, lon, var)

        # Subset contains everything requested. No need to interpolate.
        if np.all([d in dims['time'] for d in doy]) & \
                np.all([z in dims['depth'] for z in depth]) & \
                np.all([y in dims['lat'] for y in lat]) & \
                np.all([x in dims['lon'] for x in lon]):
                    dn = np.nonzero([d in doy for d in dims['time']])[0]
                    zn = np.nonzero([z in depth for z in dims['depth']])[0]
                    yn = np.nonzero([y in lat for y in dims['lat']])[0]
                    xn = np.nonzero([x in lon for x in dims['lon']])[0]
                    output = {}
                    for v in subset:
                        # output[v] = subset[v][dn, zn, yn, xn]
                        # Seriously that this is the way to do it?!!??
                        output[v] = subset[v][:, :, :, xn][:, :, yn][:, zn][dn]
                    return output

        output = {}
        for v in var:
            output[v] = ma.masked_all(
                    (doy.size, depth.size, lat.size, lon.size),
                    dtype=subset[v].dtype)

            # These interpolators don't understand Masked Arrays, but do NaN
            if subset[v].dtype in ['int32']:
                subset[v] = subset[v].astype('f')
            subset[v][ma.getmaskarray(subset[v])] = np.nan
            subset[v] = subset[v].data

        # First linear interpolate on time.
        if not (doy == dims['time']).all():
            for v in subset.keys():
                f = interp1d(dims['time'], subset[v], axis=0)
                subset[v] = f(doy)
            dims['time'] = np.atleast_1d(doy)

        if not (np.all(lat == dims['lat']) and np.all(lon == dims['lon'])):
            # Lat x Lon target coordinates are the same for all time and depth.
            points_out = []
            for latn in lat:
                for lonn in lon:
                    points_out.append([latn, lonn])
            points_out = np.array(points_out)

            # Interpolate on X/Y plane
            for v in subset:
                tmp = np.nan * np.ones(
                        (doy.size, dims['depth'].size, lat.size, lon.size),
                        dtype=subset[v].dtype)
                for nt in range(doy.size):
                    for nz in range(dims['depth'].size):
                        data = subset[v][nt, nz]
                        # The valid data
                        idx = np.nonzero(~np.isnan(data))
                        if idx[0].size > 0:
                            points = np.array([
                                dims['lat'][idx[0]], dims['lon'][idx[1]]]).T
                            values = data[idx]
                            # Interpolate along the dimensions that have more than
                            #   one position, otherwise it means that the output
                            #   is exactly on that coordinate.
                            #ind = np.array([np.unique(points[:, i]).size > 1
                            #    for i in range(points.shape[1])])
                            #assert ind.any()

                            try:
                                values_out = griddata(
                                    #np.atleast_1d(np.squeeze(points[:, ind])),
                                    np.atleast_1d(np.squeeze(points)),
                                    values,
                                    #np.atleast_1d(np.squeeze(points_out[:, ind])))
                                    np.atleast_1d(np.squeeze(points_out)))
                            except:
                                values_out = []
                                for p in points_out:
                                    try:
                                        values_out.append(griddata(
                                            np.atleast_1d(np.squeeze(points)),
                                            values,
                                            np.atleast_1d(np.squeeze(
                                                p))))
                                    except:
                                        values_out.append(np.nan)
                                values_out = np.array(values_out)

                            # Remap the interpolated value back into a 4D array
                            idx = np.isfinite(values_out)
                            for [y, x], out in zip(
                                    points_out[idx], values_out[idx]):
                                tmp[nt, nz, y==lat, x==lon] = out
                subset[v] = tmp

        # Interpolate on z
        same_depth = (np.shape(depth) == dims['depth'].shape) and \
                np.allclose(depth, dims['depth'])
        if not same_depth:
            for v in list(subset.keys()):
                try:
                    f = interp1d(dims['depth'], subset[v], axis=1, bounds_error=False)
                    # interp1d does not handle Masked Arrays
                    subset[v] = f(np.array(depth))
                except:
                    print("Fail to interpolate '%s' in depth" % v)
                    del(subset[v])

        for v in subset:
            if output[v].dtype in ['int32']:
                subset[v] = np.round(subset[v])
            output[v][:] = ma.fix_invalid(subset[v][:])

        return output
172685
def interpolation_bilinear(x, y, x1, x2, y1, y2, z11, z21, z22, z12):
    '''
    
    '''
    t = (x - x1) / (x2 - x1)
    s = (y - y1) / (y2 - y1)
    v1 = (1.0 - t) * (1.0 - s) * z11
    v2 = t * (1.0 - s) * z21
    v3 = t * s * z22
    v4 = (1.0 - t) * s * z12
    ret = v1 + v2 + v3 + v4
    return ret
174016
def interpolate_jagged(xyz,nseg):
    
    
    # Spherical coordinates specifying the angles of all line
    # segments that make up the section path
    (r,theta,phi) = sequential_spherical(xyz)
    
    # cumulative length of section path at each coordinate
    rcum = np.append(0,np.cumsum(r))

    # breakpoints for segment paths along section path
    breakpoints = np.linspace(0,rcum[-1],nseg+1)
    np.delete(breakpoints,0)
    
    # Find segment paths
    seg_paths = []
    for a in range(nseg):
        path = []
        
        # find (x,y,z) starting coordinate of path
        if a == 0:
            start_coord = xyz[0,:]
        else:
            start_coord = end_coord # start at end of last path
        path.append(start_coord)

        # find all coordinates between the start and end points
        start_length = breakpoints[a]
        end_length = breakpoints[a+1]
        mid_boolean = (rcum > start_length) & (rcum < end_length)
        mid_indices = np.nonzero(mid_boolean)[0]
        for mi in mid_indices:
            path.append(xyz[mi,:])

        # find (x,y,z) ending coordinate of path
        end_coord = find_coord(end_length,xyz,rcum,theta,phi)
        path.append(end_coord)

        # Append path to list of segment paths
        seg_paths.append(np.array(path))
    
    # Return all segment paths
    return seg_paths
1014272
def __generate_spline(self, x, y, closed=False, steps=20):
        

        if closed:
            x = x.tolist()
            x.insert(0, x[-1])
            x.append(x[1])
            x.append(x[2])

            y = y.tolist()
            y.insert(0, y[-1])
            y.append(y[1])
            y.append(y[2])

        points = np.vstack((x,y)).T

        curve = []

        if not closed:
            curve.append(points[0])

        for j in range(1, len(points)-2):
            for s in range(steps):
                t = 1. * s / steps
                p0, p1, p2, p3 = points[j-1], points[j], points[j+1], points[j+2]
                pnew = 0.5 *((2 * p1) + (-p0 + p2) * t + (2*p0 - 5*p1 + 4*p2 - p3) * t**2 + (-p0 + 3*p1- 3*p2 + p3) * t**3)
                curve.append(pnew)

        if not closed:
            curve.append(points[-1])

        curve = np.array(curve)
        return curve[:, 0], curve[:, 1]
266722
function (latlngs) {
		var result = [],
		    flat = L.Polyline._flat(latlngs);

		for (var i = 0, len = latlngs.length; i < len; i++) {
			if (flat) {
				result[i] = L.latLng(latlngs[i]);
				this._bounds.extend(result[i]);
			} else {
				result[i] = this._convertLatLngs(latlngs[i]);
			}
		}

		return result;
	}
762642
void interpolateDeformedPoint(float v_x , float v_y , Point2D_F32 deformed ) {

		// sample the closest point and x+1,y+1
		int x0 = (int)v_x;
		int y0 = (int)v_y;
		int x1 = x0+1;
		int y1 = y0+1;

		// make sure the 4 sample points are in bounds
		if( x1 >= gridCols )
			x1 = gridCols-1;
		if( y1 >= gridRows )
			y1 = gridRows-1;

		// weight along each axis
		float ax = v_x - x0;
		float ay = v_y - y0;

		// bilinear weight for each sample point
		float w00 = (1.0f - ax) * (1.0f - ay);
		float w01 = ax * (1.0f - ay);
		float w11 = ax * ay;
		float w10 = (1.0f - ax) * ay;

		// apply weights to each sample point
		Point2D_F32 d00 = getGrid(y0,x0).deformed;
		Point2D_F32 d01 = getGrid(y0,x1).deformed;
		Point2D_F32 d10 = getGrid(y1,x0).deformed;
		Point2D_F32 d11 = getGrid(y1,x1).deformed;

		deformed.set(0,0);
		deformed.x += w00 * d00.x;
		deformed.x += w01 * d01.x;
		deformed.x += w11 * d11.x;
		deformed.x += w10 * d10.x;

		deformed.y += w00 * d00.y;
		deformed.y += w01 * d01.y;
		deformed.y += w11 * d11.y;
		deformed.y += w10 * d10.y;
	}