/*
* Copyright 1998-2014 University Corporation for Atmospheric Research/Unidata
*
* Portions of this software were developed by the Unidata Program at the
* University Corporation for Atmospheric Research.
*
* Access and use of this software shall impose the following obligations
* and understandings on the user. The user is granted the right, without
* any fee or cost, to use, copy, modify, alter, enhance and distribute
* this software, and any derivative works thereof, and its supporting
* documentation for any purpose whatsoever, provided that this entire
* notice appears in all copies of the software, derivative works and
* supporting documentation. Further, UCAR requests that the user credit
* UCAR/Unidata in any publications that result from the use of this
* software or in any product that includes this software. The names UCAR
* and/or Unidata, however, may not be used in any advertising or publicity
* to endorse or promote any products or commercial entity unless specific
* written permission is obtained from UCAR/Unidata. The user also
* understands that UCAR/Unidata is not obligated to provide the user with
* any support, consulting, training or assistance of any kind with regard
* to the use, operation and performance of this software nor to provide
* the user with any updates, revisions, new versions or "bug fixes."
*
* THIS SOFTWARE IS PROVIDED BY UCAR/UNIDATA "AS IS" AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL UCAR/UNIDATA BE LIABLE FOR ANY SPECIAL,
* INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING
* FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
* NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
* WITH THE ACCESS, USE OR PERFORMANCE OF THIS SOFTWARE.
*/
package ucar.nc2.dataset.conv;
import ucar.ma2.*;
import ucar.nc2.Attribute;
import ucar.nc2.Dimension;
import ucar.nc2.NetcdfFile;
import ucar.nc2.Variable;
import ucar.nc2.constants.AxisType;
import ucar.nc2.constants.CDM;
import ucar.nc2.constants._Coordinate;
import ucar.nc2.dataset.*;
import ucar.nc2.util.CancelTask;
import java.io.IOException;
/**
* Cosmic data - version 1.
* Add time coordinate from global atts start_time, stop_time, assuming its linear along the vertical dimension.
*
* @author caron
* @since Jul 29, 2009
*/
public class Cosmic1Convention extends CoordSysBuilder {
/**
* @param ncfile the NetcdfFile to test
* @return true if we think this is a Zebra file.
*/
public static boolean isMine(NetcdfFile ncfile) {
// :start_time = 9.17028312E8; // double
// :stop_time = 9.170284104681826E8; // double
if ((null == ncfile.findDimension("MSL_alt"))
&& (null == ncfile.findDimension("time"))) {
return false;
}
// if (null == ncfile.findGlobalAttribute( "start_time")) return false;
// if (null == ncfile.findGlobalAttribute( "stop_time")) return false;
String center = ncfile.findAttValueIgnoreCase(null, "center", null);
return (center != null) && center.equals("UCAR/CDAAC");
}
/**
* _more_
*/
public Cosmic1Convention() {
this.conventionName = "Cosmic1";
}
/**
* _more_
*
* @param ds _more_
* @param cancelTask _more_
* @throws IOException _more_
*/
public void augmentDataset(NetcdfDataset ds,
CancelTask cancelTask) throws IOException {
Attribute leoAtt = ds.findGlobalAttribute("leoId");
if (leoAtt == null) {
if (ds.findVariable("time") == null) {
// create a time variable - assume its linear along the vertical dimension
double start = ds.readAttributeDouble(null, "start_time",
Double.NaN);
double stop = ds.readAttributeDouble(null, "stop_time",
Double.NaN);
if (Double.isNaN(start) && Double.isNaN(stop)) {
double top = ds.readAttributeDouble(null, "toptime",
Double.NaN);
double bot = ds.readAttributeDouble(null, "bottime",
Double.NaN);
this.conventionName = "Cosmic2";
if (top > bot) {
stop = top;
start = bot;
} else {
stop = bot;
start = top;
}
}
Dimension dim = ds.findDimension("MSL_alt");
Variable dimV = ds.findVariable("MSL_alt");
Array dimU = dimV.read();
int inscr = (dimU.getFloat(1) - dimU.getFloat(0)) > 0
? 1
: 0;
int n = dim.getLength();
double incr = (stop - start) / n;
String timeUnits = "seconds since 1980-01-06 00:00:00";
Variable timeVar = new VariableDS(ds, null, null, "time",
DataType.DOUBLE, dim.getShortName(),
timeUnits, null);
ds.addVariable(null, timeVar);
timeVar.addAttribute(new Attribute(CDM.UNITS, timeUnits));
timeVar.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Time.toString()));
int dir = ds.readAttributeInteger(null, "irs", 1);
ArrayDouble.D1 data =
(ArrayDouble.D1) Array.factory(DataType.DOUBLE,
new int[]{n});
if (inscr == 0) {
if (dir == 1) {
for (int i = 0; i < n; i++) {
data.set(i, start + i * incr);
}
} else {
for (int i = 0; i < n; i++) {
data.set(i, stop - i * incr);
}
}
} else {
for (int i = 0; i < n; i++) {
data.set(i, stop - i * incr);
}
}
timeVar.setCachedData(data, false);
}
Variable v = ds.findVariable("Lat");
if (v == null) {
v = ds.findVariable("GEO_lat");
}
v.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Lat.toString()));
Variable v1 = ds.findVariable("Lon");
if (v1 == null) {
v1 = ds.findVariable("GEO_lon");
}
v1.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Lon.toString()));
} else {
Dimension dim = ds.findDimension("time");
int n = dim.getLength();
Variable latVar = new VariableDS(ds, null, null, "Lat",
DataType.FLOAT, dim.getShortName(),
"degree", null);
latVar.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Lat.toString()));
ds.addVariable(null, latVar);
Variable lonVar = new VariableDS(ds, null, null, "Lon",
DataType.FLOAT, dim.getShortName(),
"degree", null);
lonVar.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Lon.toString()));
ds.addVariable(null, lonVar);
Variable altVar = new VariableDS(ds, null, null, "MSL_alt",
DataType.FLOAT, dim.getShortName(),
"meter", null);
altVar.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Height.toString()));
ds.addVariable(null, altVar);
// cal data array
ArrayFloat.D1 latData =
(ArrayFloat.D1) Array.factory(DataType.FLOAT,
new int[]{n});
ArrayFloat.D1 lonData =
(ArrayFloat.D1) Array.factory(DataType.FLOAT,
new int[]{n});
ArrayFloat.D1 altData =
(ArrayFloat.D1) Array.factory(DataType.FLOAT,
new int[]{n});
ArrayDouble.D1 timeData =
(ArrayDouble.D1) Array.factory(DataType.DOUBLE,
new int[]{n});
this.conventionName = "Cosmic3";
int iyr = ds.readAttributeInteger(null, "year", 2009);
int mon = ds.readAttributeInteger(null, "month", 0);
int iday = ds.readAttributeInteger(null, "day", 0);
int ihr = ds.readAttributeInteger(null, "hour", 0);
int min = ds.readAttributeInteger(null, "minute", 0);
int sec = ds.readAttributeInteger(null, "second", 0);
double start = ds.readAttributeDouble(null, "startTime",
Double.NaN);
double stop = ds.readAttributeDouble(null, "stopTime",
Double.NaN);
double incr = (stop - start) / n;
int t = 0;
// double julian = juday(mon, iday, iyr);
// cal the dtheta based pm attributes
double dtheta = gast(iyr, mon, iday, ihr, min, sec, t);
Variable tVar = ds.findVariable("time");
String timeUnits = "seconds since 1980-01-06 00:00:00"; //dtime.getUnit().toString();
tVar.removeAttributeIgnoreCase(CDM.VALID_RANGE);
tVar.removeAttributeIgnoreCase(CDM.UNITS);
tVar.addAttribute(new Attribute(CDM.UNITS, timeUnits));
tVar.addAttribute(new Attribute(_Coordinate.AxisType,
AxisType.Time.toString()));
Variable v = ds.findVariable("xLeo");
Array xLeo = v.read();
v = ds.findVariable("yLeo");
Array yLeo = v.read();
v = ds.findVariable("zLeo");
Array zLeo = v.read();
double a = 6378.1370;
double b = 6356.7523142;
IndexIterator iiter0 = xLeo.getIndexIterator();
IndexIterator iiter1 = yLeo.getIndexIterator();
IndexIterator iiter2 = zLeo.getIndexIterator();
int i = 0;
while (iiter0.hasNext()) {
double[] v_inertial = new double[3];
v_inertial[0] = iiter0.getDoubleNext(); //.getDouble(i); //.nextDouble();
v_inertial[1] = iiter1.getDoubleNext(); //.getDouble(i); //.nextDouble();
v_inertial[2] = iiter2.getDoubleNext(); //.getDouble(i); //.nextDouble();
double[] uvz = new double[3];
uvz[0] = 0.0;
uvz[1] = 0.0;
uvz[2] = 1.0;
// v_ecef should be in the (approximate) ECEF frame
// double[] v_ecf = execute(v_inertial, julian);
double[] v_ecf = spin(v_inertial, uvz, -1 * dtheta);
// cal lat/lon here
// double [] llh = ECFtoLLA(v_ecf[0]*1000, v_ecf[1]*1000, v_ecf[2]*1000, a, b);
double[] llh = xyzell(a, b, v_ecf);
latData.set(i, (float) llh[0]);
lonData.set(i, (float) llh[1]);
altData.set(i, (float) llh[2]);
timeData.set(i, start + i * incr);
i++;
}
latVar.setCachedData(latData, false);
lonVar.setCachedData(lonData, false);
altVar.setCachedData(altData, false);
tVar.setCachedData(timeData, false);
}
ds.finish();
}
/* private DateTime getDateTime(int year, int month, int day, int hour,
int min, int sec) throws Exception {
GregorianCalendar convertCal =
new GregorianCalendar(DateUtil.TIMEZONE_GMT);
convertCal.clear();
convertCal.set(Calendar.YEAR, year);
//The MONTH is 0 based. The incoming month is 1 based
convertCal.set(Calendar.MONTH, month - 1);
convertCal.set(Calendar.DAY_OF_MONTH, day);
convertCal.set(Calendar.HOUR_OF_DAY, hour);
convertCal.set(Calendar.MINUTE, min);
convertCal.set(Calendar.SECOND, sec);
return new DateTime(convertCal.getTime());
} */
/**
* _more_
*
* @param ncDataset _more_
* @param v _more_
* @return _more_
*/
protected AxisType getAxisType(NetcdfDataset ncDataset,
VariableEnhanced v) {
String name = v.getShortName();
if (name.equals("time")) {
return AxisType.Time;
}
if (name.equals("Lat") || name.equals("GEO_lat")) {
return AxisType.Lat;
}
if (name.equals("Lon") || name.equals("GEO_lon")) {
return AxisType.Lon;
}
// if (name.equals("xLeo") ) return AxisType.GeoX;
// if (name.equals("yLeo") ) return AxisType.GeoY;
if (name.equals("MSL_alt")) {
return AxisType.Height;
}
return null;
}
/**
* NAME : XYZELL
* <p/>
* CALL XYZELL(A,B,XSTAT,XSTELL)
* <p/>
* PURPOSE : COMPUTATION OF ELLIPSOIDAL COORDINATES "XSTELL"
* GIVEN THE CARTESIAN COORDINATES "XSTAT"
* <p/>
* PARAMETERS :
* IN : A : SEMI-MAJOR AXIS OF THE REFERENCE R*8
* ELLIPSOID IN METERS
* B : SEMI-MINOR AXIS OF THE REFERENCE R*8
* ELLIPSOID IN METERS
* DXELL(3): TRANSLATION COMPONENTS FROM THE R*8
* ORIGIN OF THE CART. COORD. SYSTEM
* (X,Y,Z) TO THE CENTER OF THE REF.
* ELLIPSOID (IN METRES)
* SCELL : SCALE FACTOR BETWEEN REF. ELLIPSOID R*8
* AND WGS-84
* XSTAT(3): CARTESIAN COORDINATES (M) R*8
* OUT : XSTELL(3): ELLIPSOIDAL COORDINATES R*8
* XSTELL(1): ELL. LATITUDE (RADIAN)
* XSTELL(2): ELL. LONGITUDE (RADIAN)
* XSTELL(3): ELL. HEIGHT (M)
* <p/>
* SR CALLED : DMLMTV
* <p/>
* REMARKS : ---
* <p/>
* AUTHOR : M. ROTHACHER
* <p/>
* VERSION : 3.4 (JAN 93)
* <p/>
* CREATED : 87/11/03 12:32 LAST MODIFIED : 88/11/21 17:36
* <p/>
* COPYRIGHT : ASTRONOMICAL INSTITUTE
* 1987 UNIVERSITY OF BERNE
* SWITZERLAND
*
* @param a _more_
* @param b _more_
* @param xstat _more_
* @return _more_
*/
public double[] xyzell(double a, double b, double[] xstat) {
double[] xstell = new double[3];
double e2, s, rlam, zps, h, phi, n, hp, phip;
int i, niter;
e2 = (a * a - b * b) / (a * a);
s = Math.sqrt(xstat[0] * xstat[0] + xstat[1] * xstat[1]);
rlam = Math.atan2(xstat[1], xstat[0]);
zps = xstat[2] / s;
h = Math.sqrt(xstat[0] * xstat[0] + xstat[1] * xstat[1]
+ xstat[2] * xstat[2]) - a;
phi = Math.atan(zps / (1.0 - e2 * a / (a + h)));
niter = 0;
for (i = 1; i <= 10000000; i++) {
n = a / Math.sqrt(1.0 - e2 * Math.sin(phi) * Math.sin(phi));
hp = h;
phip = phi;
h = s / Math.cos(phi) - n;
phi = Math.atan(zps / (1.0 - e2 * n / (n + h)));
niter = niter + 1;
if ((Math.abs(phip - phi) <= 1.e-11)
&& (Math.abs(hp - h) <= 1.e-5)) {
break;
}
if (niter >= 10) {
phi = -999.0;
rlam = -999.0;
h = -999.0;
break;
}
}
xstell[0] = phi * 180 / 3.1415926;
xstell[1] = rlam * 180 / 3.1415926;
xstell[2] = h;
return xstell;
}
/**
* gast.f
* Compute hour angle dtheta
*
* ! iyr, mon, iday, ihr, min and sec form a base (epoch) time,
* ! t is an offset from the base time in seconds
* ! dtheta is the output hour angle in radians
*
* Calculation of local time
*
* ! glon -- East longitude in degrees, -180 to 180
*
*
* call vprod.f spin.f rnorm.f
* Calculation of the unit vector normal to the occultation plane
* (clockwise rotated from GPS to LEO)
*
* @param iyr _more_
* @param imon _more_
* @param iday _more_
* @param ihr _more_
* @param imin _more_
* @param sec _more_
* @param dsec _more_
*
* @return _more_
*/
/*
double dtheta = gast(iyr,mon,iday,ihr,min,sec,t)
utc=ihr*1.d0+min/60.d0
timloc=utc+24.d0*glon/360.d0
if (timloc.gt.24.d0) timloc=timloc-24.d0
if (timloc.lt.0.d0) timloc=timloc+24.d0
*/
// In the inertial reference frame
/*
v_inertial(1)= ! Inertial GPS position vectors, XYZ
v_inertial(2)=
v_inertial(3)=
*/
// In the Earth-fixed reference frame
// Z axis to rotate around (unit vector Z)
/*
uvz(1)=0.0;
uvz(2)=0.0;
uvz(3)=1.0;
double [] v_ecf = spin(v_inertial,uvz,-180.d0*dtheta/pi)
*/
// after this call, v_ecef should be in the (approximate) ECEF frame
/**
* ----------------------------------------------------------------------
* gast.f
* <p/>
* This subroutine computes the Greenwich Apparent Siderial
* Time angle given a UTC date and time.
* <p/>
* parameter Input parameters:
* Inputs:
*
* @param iyr integer year 1995
* @param imon integer month 5
* @param iday integer day 5
* @param ihr integer hour 5
* @param imin integer minute 5
* @param sec double second 31.0
* @param dsec double second 0.0
* Outputs:
* @return theta GAST angle in radians
* <p/>
* author Bill Schreiner
* @since May 1995
* version $URL: svn://ursa.cosmic.ucar.edu/trunk/src/roam/gast.f $ $Id: gast.f 10129 2008-07-30 17:10:52Z dhunt $
* -----------------------------------------------------------------------
*/
public double gast(int iyr, int imon, int iday, int ihr, int imin,
double sec, double dsec) {
//
// implicit double precision (a-h,o-z)
// character(len=*), parameter :: header = '$URL: svn://ursa.cosmic.ucar.edu/trunk/src/roam/gast.f $ $Id: gast.f 10129 2008-07-30 17:10:52Z dhunt $'
//
// Coordinate transform from the celestial inertial reference frame to the geo-
// centered Greenwich reference frame.
// Call a subroutine to calculate the Julian day "djd":
double djd = juday(imon, iday, iyr); //djd=julean day.
double tu = (djd - 2451545.0) / 36525.0;
double gmst = 24110.548410 + 8640184.8128660 * tu
+ 0.093104 * tu * tu - 6.2E-6 * Math.pow(tu, 3); // !gmst=Greenwich mean...
double utco = (ihr * 3600) + (imin * 60) + sec;
return togreenw(dsec, utco, gmst);
}
/**
* JDAY calculates the Julian Day number (JD) from the Gregorian month
* ,day, and year (M,D,Y). (NOT VALID BEFORE 10/15/1582)
*
* @param M _more_
* @param D _more_
* @param Y _more_
* @return _more_
*/
public double juday(int M, int D, int Y) {
double JD;
double IY = Y - (12 - M) / 10;
double IM = M + 1 + 12 * ((12 - M) / 10);
double I = IY / 100;
double J = 2 - I + I / 4 + Math.round(365.25 * IY)
+ Math.round(30.6001 * IM);
JD = (J + D + 1720994.50);
return JD;
}
/**
* This subroutine is to transform the locations and velocities of the GPS and
* LEO satellites from the celestial inertial reference frame to the Earth
* centered Greenwich reference frame.
* The dummy arguments iyear, month and iday are the calender year, month and
* day of the occultation event. The dummy arguments ihour, minute and sec
* are the UTC time.
* Reference: Astronomical Alamanus, 1993
* <p/>
* Modified subroutine from Dasheng's code.
*
* @param rectt _more_
* @param utco _more_
* @param gmst _more_
* @return _more_
*/
public double togreenw(double rectt, double utco, double gmst) {
double pi = Math.acos(-1.00);
//
// For each occultation ID, its TU and GMST are the same. However, every
// occultation event takes place at gmst+uts, uts is progressively increasing
// with every occultation event.
double utc = (utco + rectt) * 1.0027379093;
gmst = gmst + utc; //in seconds, without eoe correction.
// gmst may be a positive number or may be a negative number.
while (gmst < 0.0) {
gmst = gmst + 86400.00;
}
while (gmst > 86400.00) {
gmst = gmst - 86400.00;
}
// gmst = the Greenwich mean sidereal time.
// This gmst is without the corrections from the equation of equinoxes. For
// GPS/MET applications, the corrections from equation of equinoxes is not
// necessary because of the accurary needed.
return gmst * 2.0 * pi / 86400.0; //!*** This is the THETA in radian.
}
/**
* ----------------------------------------------------------------------
* file spin.f
* <p/>
* This subroutine rotates vector V1 around vector VS
* at angle A. V2 is the vector after the rotation.
* <p/>
* <p/>
* parameter Input parameters:
* v1 - Vector to be rotated
* vs - Vector around which to rotate v1
* a - angle of rotation
* Output parameters:
* v2 - output vector
* <p/>
* S.V.Sokolovskiy
* URL: svn://ursa.cosmic.ucar.edu/trunk/src/roam/spin.f $ $Id: spin.f 10129 2008-07-30 17:10:52Z dhunt $
* -----------------------------------------------------------------------
*
* @param v1 - Vector to be rotated
* @param vs - Vector around which to rotate v1
* @param a - angle of rotation
* @return _more_
*/
public double[] spin(double[] v1, double[] vs, double a) {
// implicit real*8(a-h,o-z)
// dimension v1(3),vs(3),vsn(3),v2(3),v3(3),s(3,3)
// Calculation of the unit vector around which
// the rotation should be done.
double[] v2 = new double[3];
double[] vsn = new double[3];
double[] v3 = new double[3];
double vsabs = Math.sqrt(vs[0] * vs[0] + vs[1] * vs[1]
+ vs[2] * vs[2]);
for (int i = 0; i < 3; i++) {
vsn[i] = vs[i] / vsabs;
}
// Calculation of the rotation matrix.
double a1 = Math.cos(a);
double a2 = 1.0 - a1;
double a3 = Math.sin(a);
double[][] s = new double[3][3];
s[0][0] = a2 * vsn[0] * vsn[0] + a1;
s[0][1] = a2 * vsn[0] * vsn[1] - a3 * vsn[2];
s[0][2] = a2 * vsn[0] * vsn[2] + a3 * vsn[1];
s[1][0] = a2 * vsn[1] * vsn[0] + a3 * vsn[2];
s[1][1] = a2 * vsn[1] * vsn[1] + a1;
s[1][2] = a2 * vsn[1] * vsn[2] - a3 * vsn[0];
s[2][0] = a2 * vsn[2] * vsn[0] - a3 * vsn[1];
s[2][1] = a2 * vsn[2] * vsn[1] + a3 * vsn[0];
s[2][2] = a2 * vsn[2] * vsn[2] + a1;
// Calculation of the rotated vector.
for (int i = 0; i < 3; i++) {
v3[i] = s[i][0] * v1[0] + s[i][1] * v1[1] + s[i][2] * v1[2];
}
System.arraycopy(v3, 0, v2, 0, 3);
return v2;
}
/**
* _more_
*/
protected final static double RTD = 180. / Math.PI;
/**
* _more_
*/
protected final static double DTR = Math.PI / 180.;
/**
* _more_
*
* @param eci _more_
* @param julian _more_
* @return _more_
*/
public double[] execute(double[] eci, double julian) {
double Xi = eci[0];
double Yi = eci[1];
double Zi = eci[2];
double c, s;
double GHA;
double[] ecef = new double[3];
//Compute GHAD
/* System generated locals */
double d__1, d__2, d__3;
/* Local variables */
double tsec, tday, gmst, t, omega, tfrac, tu, dat;
/* INPUT IS TIME "secondsSince1970" IN SECONDS AND "TDAY" */
/* WHICH IS WHOLE DAYS FROM 1970 JAN 1 0H */
/* THE OUTPUT IS GREENWICH HOUR ANGLE IN DEGREES */
/* XOMEGA IS ROTATION RATE IN DEGREES/SEC */
/* FOR COMPATABILITY */
tday = (double) ((int) (julian / 86400.));
tsec = julian - tday * 86400;
/* THE NUMBER OF DAYS FROM THE J2000 EPOCH */
/* TO 1970 JAN 1 0H UT1 IS -10957.5 */
t = tday - (float) 10957.5;
tfrac = tsec / 86400.;
dat = t;
tu = dat / 36525.;
/* Computing 2nd power */
d__1 = tu;
/* Computing 3rd power */
d__2 = tu;
d__3 = d__2;
gmst = tu * 8640184.812866 + 24110.54841 + d__1 * d__1 * .093104
- d__3 * (d__2 * d__2) * 6.2e-6;
/* COMPUTE THE EARTH'S ROTATION RATE */
/* Computing 2nd power */
d__1 = tu;
omega = tu * 5.098097e-6 + 86636.55536790872 - d__1 * d__1 * 5.09e-10;
/* COMPUTE THE GMST AND GHA */
// da is earth nutation - currently unused
double da = 0.0;
gmst = gmst + omega * tfrac + da * RTD * 86400. / 360.;
gmst = gmst % 86400;
if (gmst < 0.) {
gmst += 86400.;
}
gmst = gmst / 86400. * 360.;
//ghan = gmst;
// returns gha in radians
gmst = gmst * DTR;
GHA = gmst;
//RotateZ
c = Math.cos(GHA);
s = Math.sin(GHA);
double X = c * Xi + s * Yi;
double Y = -s * Xi + c * Yi;
//Set outputs
ecef[0] = X;
ecef[1] = Y;
ecef[2] = Zi;
return ecef;
}
/**
* comparing api to others
*
* @param x _more_
* @param y _more_
* @param z _more_
* @param a _more_
* @param b _more_
* @return _more_
*/
public static double[] ECFtoLLA(double x, double y, double z,
double a, double b) {
double longitude = Math.atan2(y, x);
double ePrimeSquared = (a * a - b * b) / (b * b);
double p = Math.sqrt(x * x + y * y);
double theta = Math.atan((z * a) / (p * b));
double sineTheta = Math.sin(theta);
double cosTheta = Math.cos(theta);
double f = 1 / 298.257223563;
double e2 = 2 * f - f * f;
double top = z + ePrimeSquared * b * sineTheta * sineTheta
* sineTheta;
double bottom = p - e2 * a * cosTheta * cosTheta * cosTheta;
double geodeticLat = Math.atan(top / bottom);
double sineLat = Math.sin(geodeticLat);
double N = a / Math.sqrt(1 - e2 * sineLat * sineLat);
double altitude = (p / Math.cos(geodeticLat)) - N;
// maintain longitude btw -PI and PI
if (longitude > Math.PI) {
longitude -= 2 * Math.PI;
} else if (longitude < -Math.PI) {
longitude += 2 * Math.PI;
}
return new double[]{geodeticLat, longitude, altitude};
}
}