List of usage examples for org.apache.commons.math.util FastMath sin
public static double sin(double x)
From source file:cz.paulrz.montecarlo.random.Sobol.java
@Override public double nextGaussian() { final double random; if (Double.isNaN(nextGaussian)) { // generate a new pair of gaussian numbers final double[] xs = nextPoint(); final double x = xs[0]; final double y = xs[1]; final double alpha = 2 * FastMath.PI * x; final double r = FastMath.sqrt(-2 * FastMath.log(y)); random = r * FastMath.cos(alpha); nextGaussian = r * FastMath.sin(alpha); } else {//from w w w.jav a 2 s .c o m // use the second element of the pair already generated random = nextGaussian; nextGaussian = Double.NaN; } return random; }
From source file:edu.berkeley.path.bots.core.Coordinate.java
public double distanceHaversineInMeters(Coordinate otherCoord) { // SRID's must match and not be null. // This can't be null, duh, and we don't check if otherCoord is, which // means the NullPointerException will be thrown. if (null == this.srid() || null == otherCoord.srid()) { throw new ClassCastException("This distance function uses the spheroid distance, but you " + "are using Coordinates with null SRID~s (from a PostgreSQL " + "point?). This doesn't really make any sense and you " + "probably want to use .distanceCarteasianInSRUnits( other )" + "instead."); } else if (!this.srid().equals(otherCoord.srid())) { throw new ClassCastException("The SRID of otherCoord does't match this one."); }/*from w ww . j av a 2 s . c o m*/ final double piOver180 = Math.PI / 180.0; final double earthRadiusInMeters = 6367000; final double lon1 = piOver180 * this.lon(); final double lat1 = piOver180 * this.lat(); final double lon2 = piOver180 * otherCoord.lon(); final double lat2 = piOver180 * otherCoord.lat(); // Haversine formula: final double dlon = lon2 - lon1; final double dlat = lat2 - lat1; final double a = FastMath.sin(dlat / 2) * FastMath.sin(dlat / 2) + FastMath.cos(lat1) * FastMath.cos(lat2) * FastMath.sin(dlon / 2) * FastMath.sin(dlon / 2); final double c = 2 * FastMath.atan2(FastMath.sqrt(a), FastMath.sqrt(1 - a)); return earthRadiusInMeters * c; }
From source file:edu.berkeley.path.bots.core.Coordinate.java
/** * Implementation of Thaddeus Vincenty's algorithms to solve the direct and * inverse geodetic problems.//w w w . ja va 2s . co m * <p/> * This implementation may be as good as the PostGIS provided one, the * errors seem to be pretty small (10 ^ -3 meters), but a thorough * check/analysis has not been done. * <p/> * For more information, see Vincenty's original publication on the NOAA * web-site: <a href="http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf"> * http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf</a> * * @param otherCoord * @return the distance calculated this way. * @throws ClassCastException * if the SRID~s don't match or are null. * @throws NullPointerException * if otherCoord is null. */ public double distanceVincentyInMeters(Coordinate otherCoord) { // SRID's must match and not be null. // This can't be null, duh, and we don't check if otherCoord is, which // means the NullPointerException will be thrown. if (null == this.srid() || null == otherCoord.srid()) { throw new ClassCastException("This distance function uses the spheroid distance, but you " + "are using Coordinates with null SRID~s (from a PostgreSQL " + "point?). This doesn't really make any sense and you " + "probably want to use .distanceCarteasianInSRUnits( other )" + "instead."); } else if (!this.srid().equals(otherCoord.srid())) { throw new ClassCastException("The SRID of otherCoord does't match this one."); } // else we are sure they are the same, and not null. final double a; final double b; final double f; // Set a, b, and f. switch (this.srid()) { case 4326: a = 6378137.0; b = 6356752.3142; f = (a - b) / a; break; default: throw new ClassCastException("The given SRID is not entered in this function, you can" + " get it by running the DB query defined in a function" + " in the NETCONFIG project called " + "util.SpatialDatabaseQueries.getSpheroidStringCached()" + ". This should return a spheroid string like:\n" + "\tSPHEROID(\"WGS 84\",6378137,298.257223563)\n" + "which is SPHEROID( desc, d1, d2 ), so then\n" + "\ta = d1\n" + "\tb = ( d2 * d1 - d1 ) / d2\n" + "\tf = 1 / d2\n" + "Then you can add the numbers to the switch statement" + " in the Java."); } // CHECKSTYLE:OFF // Now a, b, and f should be set. // // All constants below refer to Vincenty's publication (see above). // // get parameters as radians final double PiOver180 = Math.PI / 180.0; final double phi1 = PiOver180 * this.lat(); final double lambda1 = PiOver180 * this.lon(); final double phi2 = PiOver180 * otherCoord.lat(); final double lambda2 = PiOver180 * otherCoord.lon(); // calculations final double a2 = a * a; final double b2 = b * b; final double a2b2b2 = (a2 - b2) / b2; final double omega = lambda2 - lambda1; final double tanphi1 = FastMath.tan(phi1); final double tanU1 = (1.0 - f) * tanphi1; final double U1 = FastMath.atan(tanU1); final double sinU1 = FastMath.sin(U1); final double cosU1 = FastMath.cos(U1); final double tanphi2 = FastMath.tan(phi2); final double tanU2 = (1.0 - f) * tanphi2; final double U2 = FastMath.atan(tanU2); final double sinU2 = FastMath.sin(U2); final double cosU2 = FastMath.cos(U2); final double sinU1sinU2 = sinU1 * sinU2; final double cosU1sinU2 = cosU1 * sinU2; final double sinU1cosU2 = sinU1 * cosU2; final double cosU1cosU2 = cosU1 * cosU2; // Below are the re-assignable fields // eq. 13 double lambda = omega; // intermediates we'll need to compute 's' double A = 0.0; double B = 0.0; double sigma = 0.0; double deltasigma = 0.0; double lambda0; final int num_iters = 5; // 20 final double change_threshold = 1e-5; // 0.0000000000001; for (int i = 0; i < num_iters; i++) { lambda0 = lambda; final double sinlambda = FastMath.sin(lambda); final double coslambda = FastMath.cos(lambda); // eq. 14 final double sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) + (cosU1sinU2 - sinU1cosU2 * coslambda) * (cosU1sinU2 - sinU1cosU2 * coslambda); final double sinsigma = FastMath.sqrt(sin2sigma); // eq. 15 final double cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda); // eq. 16 sigma = FastMath.atan2(sinsigma, cossigma); // eq. 17 Careful! sin2sigma might be almost 0! final double sinalpha = (sin2sigma == 0) ? 0.0 : cosU1cosU2 * sinlambda / sinsigma; final double alpha = FastMath.asin(sinalpha); final double cosalpha = FastMath.cos(alpha); final double cos2alpha = cosalpha * cosalpha; // eq. 18 Careful! cos2alpha might be almost 0! final double cos2sigmam = cos2alpha == 0.0 ? 0.0 : cossigma - 2 * sinU1sinU2 / cos2alpha; final double u2 = cos2alpha * a2b2b2; final double cos2sigmam2 = cos2sigmam * cos2sigmam; // eq. 3 A = 1.0 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))); // eq. 4 B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2))); // eq. 6 deltasigma = B * sinsigma * (cos2sigmam + B / 4 * (cossigma * (-1 + 2 * cos2sigmam2) - B / 6 * cos2sigmam * (-3 + 4 * sin2sigma) * (-3 + 4 * cos2sigmam2))); // eq. 10 final double C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha)); // eq. 11 (modified) lambda = omega + (1 - C) * f * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam2))); // see how much improvement we got final double change = FastMath.abs((lambda - lambda0) / lambda); if ((i > 1) && (change < change_threshold)) { break; } } // for // eq. 19 return b * A * (sigma - deltasigma); // CHECKSTYLE:ON }
From source file:org.esa.beam.framework.datamodel.TiePointGrid.java
private synchronized void initDiscont() { TiePointGrid base = this; final float[] tiePoints = base.getTiePoints(); final float[] sinTiePoints = new float[tiePoints.length]; final float[] cosTiePoints = new float[tiePoints.length]; for (int i = 0; i < tiePoints.length; i++) { float tiePoint = tiePoints[i]; sinTiePoints[i] = (float) FastMath.sin(MathUtils.DTOR * tiePoint); cosTiePoints[i] = (float) FastMath.cos(MathUtils.DTOR * tiePoint); }/*from ww w .ja va 2s . c o m*/ sinGrid = new TiePointGrid(base.getName(), base.getRasterWidth(), base.getRasterHeight(), base.getOffsetX(), base.getOffsetY(), base.getSubSamplingX(), base.getSubSamplingY(), sinTiePoints); cosGrid = new TiePointGrid(base.getName(), base.getRasterWidth(), base.getRasterHeight(), base.getOffsetX(), base.getOffsetY(), base.getSubSamplingX(), base.getSubSamplingY(), cosTiePoints); }
From source file:org.esa.nest.eo.GeoUtils.java
/** * Convert geodetic coordinate into cartesian XYZ coordinate with specified geodetic system. * @param latitude The latitude of a given pixel (in degree). * @param longitude The longitude of the given pixel (in degree). * @param altitude The altitude of the given pixel (in m) * @param xyz The xyz coordinates of the given pixel. * @param geoSystem The geodetic system. *//*from ww w . j av a 2 s .co m*/ public static void geo2xyz(final double latitude, final double longitude, final double altitude, final double xyz[], final EarthModel geoSystem) { double a = 0.0; double e2 = 0.0; if (geoSystem == EarthModel.WGS84) { a = WGS84.a; e2 = WGS84.e2; } else if (geoSystem == EarthModel.GRS80) { a = GRS80.a; e2 = GRS80.e2; } else { throw new OperatorException("Incorrect geodetic system"); } final double lat = latitude * org.esa.beam.util.math.MathUtils.DTOR; final double lon = longitude * org.esa.beam.util.math.MathUtils.DTOR; final double sinLat = FastMath.sin(lat); final double cosLat = FastMath.cos(lat); final double N = a / Math.sqrt(1 - e2 * sinLat * sinLat); xyz[0] = (N + altitude) * cosLat * FastMath.cos(lon); // in m xyz[1] = (N + altitude) * cosLat * FastMath.sin(lon); // in m xyz[2] = ((1 - e2) * N + altitude) * sinLat; // in m }
From source file:org.esa.nest.eo.GeoUtils.java
/** * Convert geodetic coordinate into cartesian XYZ coordinate with specified geodetic system. * @param latitude The latitude of a given pixel (in degree). * @param longitude The longitude of the given pixel (in degree). * @param altitude The altitude of the given pixel (in m) * @param xyz The xyz coordinates of the given pixel. *//*from w w w . ja v a2s . c om*/ public static void geo2xyzWGS84(final double latitude, final double longitude, final double altitude, final double xyz[]) { final double lat = latitude * org.esa.beam.util.math.MathUtils.DTOR; final double lon = longitude * org.esa.beam.util.math.MathUtils.DTOR; final double sinLat = FastMath.sin(lat); final double N = (WGS84.a / Math.sqrt(1.0 - WGS84.e2 * sinLat * sinLat)); final double NcosLat = (N + altitude) * FastMath.cos(lat); xyz[0] = NcosLat * FastMath.cos(lon); // in m xyz[1] = NcosLat * FastMath.sin(lon); // in m xyz[2] = (N + altitude - WGS84.e2 * N) * sinLat; //xyz[2] = (WGS84.e2inv * N + altitude) * sinLat; // in m }
From source file:org.esa.nest.eo.GeoUtils.java
/** * Convert cartesian XYZ coordinate into geodetic coordinate with specified geodetic system. * @param xyz The xyz coordinate of the given pixel. * @param geoPos The geodetic coordinate of the given pixel. * @param geoSystem The geodetic system. *//*from ww w .ja v a 2 s. co m*/ public static void xyz2geo(final double xyz[], final GeoPos geoPos, final EarthModel geoSystem) { double a = 0.0; double b = 0.0; double e2 = 0.0; double ep2 = 0.0; if (geoSystem == EarthModel.WGS84) { a = WGS84.a; b = WGS84.b; e2 = WGS84.e2; ep2 = WGS84.ep2; } else if (geoSystem == EarthModel.GRS80) { a = GRS80.a; b = GRS80.b; e2 = GRS80.e2; ep2 = GRS80.ep2; } else { throw new OperatorException("Incorrect geodetic system"); } final double x = xyz[0]; final double y = xyz[1]; final double z = xyz[2]; final double s = Math.sqrt(x * x + y * y); final double theta = FastMath.atan(z * a / (s * b)); geoPos.lon = (float) (FastMath.atan(y / x) * org.esa.beam.util.math.MathUtils.RTOD); if (geoPos.lon < 0.0 && y >= 0.0) { geoPos.lon += 180.0; } else if (geoPos.lon > 0.0 && y < 0.0) { geoPos.lon -= 180.0; } geoPos.lat = (float) (FastMath .atan((z + ep2 * b * FastMath.pow(FastMath.sin(theta), 3)) / (s - e2 * a * FastMath.pow(FastMath.cos(theta), 3))) * org.esa.beam.util.math.MathUtils.RTOD); }
From source file:org.esa.nest.eo.GeoUtils.java
/** * Convert cartesian XYZ coordinate into geodetic coordinate with specified geodetic system. * @param xyz The xyz coordinate of the given pixel. * @param geoPos The geodetic coordinate of the given pixel. *//*from ww w .ja v a 2s. c o m*/ public static void xyz2geoWGS84(final double xyz[], final GeoPos geoPos) { final double x = xyz[0]; final double y = xyz[1]; final double z = xyz[2]; final double s = Math.sqrt(x * x + y * y); final double theta = FastMath.atan(z * WGS84.a / (s * WGS84.b)); geoPos.lon = (float) (FastMath.atan(y / x) * org.esa.beam.util.math.MathUtils.RTOD); if (geoPos.lon < 0.0 && y >= 0.0) { geoPos.lon += 180.0; } else if (geoPos.lon > 0.0 && y < 0.0) { geoPos.lon -= 180.0; } geoPos.lat = (float) (FastMath .atan((z + WGS84.ep2 * WGS84.b * FastMath.pow(FastMath.sin(theta), 3)) / (s - WGS84.e2 * WGS84.a * FastMath.pow(FastMath.cos(theta), 3))) * org.esa.beam.util.math.MathUtils.RTOD); }
From source file:org.esa.nest.eo.GeoUtils.java
/** * Convert polar coordinates to Cartesian vector. * <p>/*from ww w. java 2 s . com*/ * <b>Definitions:<b/> * <p>Latitude: angle from XY-plane towards +Z-axis.<p/> * <p>Longitude: angle in XY-plane measured from +X-axis towards +Y-axis.<p/> * </p> * <p> * Note: Apache's FastMath used in implementation. * </p> * @param latitude The latitude of a given pixel (in degree). * @param longitude The longitude of the given pixel (in degree). * @param radius The radius of the given point (in m) * @param xyz The return array vector of X, Y and Z coordinates for the input point. * * @author Petar Marikovic, PPO.labs */ public static void polar2cartesian(final double latitude, final double longitude, final double radius, final double xyz[]) { final double latRad = latitude * DTOR; final double lonRad = longitude * DTOR; final double sinLat = FastMath.sin(latRad); final double cosLat = FastMath.cos(latRad); xyz[0] = radius * cosLat * FastMath.cos(lonRad); xyz[1] = radius * cosLat * FastMath.sin(lonRad); xyz[2] = radius * sinLat; }
From source file:org.esa.nest.gpf.ALOSDeskewingOp.java
/** * Compute deskewing shift for pixel at (0,0). * @throws Exception The exceptions./* w w w.jav a2 s . c om*/ */ private void computeShift() throws Exception { if (useFAQShiftOnly) { return; } final stateVector v = getOrbitStateVector(firstLineTime); final double slr = slantRangeToFirstPixel + 0 * rangeSpacing; final double fd = getDopplerFrequency(0); // fractional shift final double[] lookYaw = new double[2]; computeLookYawAngles(v, slr, fd, lookYaw); fracShift = FastMath.sin(lookYaw[0]) * FastMath.sin(lookYaw[1]); if (useMapreadyShiftOnly) { return; } // compute absolute shift final String demName = "SRTM 3Sec"; final String demResamplingMethod = ResamplingFactory.BILINEAR_INTERPOLATION_NAME; DEMFactory.validateDEM(demName, sourceProduct); ElevationModel dem = DEMFactory.createElevationModel(demName, demResamplingMethod); final float demNoDataValue = dem.getDescriptor().getNoDataValue(); GeoPos geoPos = new GeoPos(); for (int y = 0; y < sourceImageHeight; y++) { sourceProduct.getGeoCoding().getGeoPos(new PixelPos(0.5f, y + 0.5f), geoPos); final double lat = geoPos.lat; double lon = geoPos.lon; if (lon >= 180.0) { lon -= 360.0; } final double alt = dem.getElevation(new GeoPos((float) lat, (float) lon)); if (alt == demNoDataValue) { continue; } final double[] earthPoint = new double[3]; final double[] sensorPos = new double[3]; GeoUtils.geo2xyzWGS84(geoPos.getLat(), geoPos.getLon(), alt, earthPoint); final double zeroDopplerTime = SARGeocoding.getEarthPointZeroDopplerTime(firstLineTime, lineTimeInterval, radarWaveLength, earthPoint, sensorPosition, sensorVelocity); if (zeroDopplerTime == SARGeocoding.NonValidZeroDopplerTime) { continue; } final double slantRange = SARGeocoding.computeSlantRange(zeroDopplerTime, timeArray, xPosArray, yPosArray, zPosArray, earthPoint, sensorPos); final double zeroDopplerTimeWithoutBias = zeroDopplerTime + slantRange / Constants.lightSpeedInMetersPerDay; absShift = (zeroDopplerTimeWithoutBias - firstLineTime) / lineTimeInterval - y; return; } absShift = computeFAQShift(v, 0); }