001 /*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2013
005 * Space Science and Engineering Center (SSEC)
006 * University of Wisconsin - Madison
007 * 1225 W. Dayton Street, Madison, WI 53706, USA
008 * https://www.ssec.wisc.edu/mcidas
009 *
010 * All Rights Reserved
011 *
012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
013 * some McIDAS-V source code is based on IDV and VisAD source code.
014 *
015 * McIDAS-V is free software; you can redistribute it and/or modify
016 * it under the terms of the GNU Lesser Public License as published by
017 * the Free Software Foundation; either version 3 of the License, or
018 * (at your option) any later version.
019 *
020 * McIDAS-V is distributed in the hope that it will be useful,
021 * but WITHOUT ANY WARRANTY; without even the implied warranty of
022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
023 * GNU Lesser Public License for more details.
024 *
025 * You should have received a copy of the GNU Lesser Public License
026 * along with this program. If not, see http://www.gnu.org/licenses.
027 */
028 package edu.wisc.ssec.mcidasv.data.adde.sgp4;
029
030 /**
031 * Earth computation functions
032 */
033 public class GeoFunctions
034 {
035
036 // same as below except the function takes Julian Date
037 /**
038 * Compute Geodetic Latatude/Longitude/Altitude from Mean of Date position vector and Date
039 *
040 * @param modPos Mean of date position vector
041 * @param mjd modified julian date (is this UTC or TT?) guessing UTC
042 * @return vector of geodetic [latitude,longitude,altitude]
043 */
044 public static double[] GeodeticLLA(double[] modPos, double mjd)
045 {
046 // calculate days since Y2k
047 // MJD = JD - 2400000.5
048 // passed in MJD - 51544.5
049 //double daysSinceY2k = (julDate - 2400000.5)-51544.5;
050 double daysSinceY2k = mjd - 51544.5;
051
052 return calculateGeodeticLLA(modPos, daysSinceY2k);
053 } // GeodeticLLA
054
055 // returned as Lat, Long, Alt
056 // LLA = corrected for time (geographical coordinates)
057 // for handling geodetic coordinates
058 // r = TEME positions, d = days since Y2K
059 private static double[] calculateGeodeticLLA(double[] r, double d)
060 {
061 double R_equ= AstroConst.R_Earth; // Equator radius [m]
062 double f = AstroConst.f_Earth; // Flattening
063
064 double eps_mach = 2.22E-16; // machine precision (double?)
065
066 final double eps = 1.0e3*eps_mach; // Convergence criterion
067 final double epsRequ = eps*R_equ;
068 final double e2 = f*(2.0-f); // Square of eccentricity
069
070 final double X = r[0]; // Cartesian coordinates
071 final double Y = r[1];
072 final double Z = r[2];
073 final double rho2 = X*X + Y*Y; // Square of distance from z-axis
074
075 // original class vars
076 // double lon;
077 // double lat;
078 //double h;
079 double[] LLA = new double[3];
080
081
082 // Check validity of input data
083 if (MathUtils.norm(r)==0.0)
084 {
085 System.out.println(" invalid input in Geodetic constructor");
086 LLA[1]=0.0;
087 LLA[0]=0.0;
088 LLA[2]=-AstroConst.R_Earth;
089 return LLA;
090 }
091
092 // Iteration
093 double dZ, dZ_new, SinPhi;
094 double ZdZ, Nh, N;
095
096 dZ = e2*Z;
097 for(;;)
098 {
099 ZdZ = Z + dZ;
100 Nh = Math.sqrt( rho2 + ZdZ*ZdZ );
101 SinPhi = ZdZ / Nh; // Sine of geodetic latitude
102 N = R_equ / Math.sqrt(1.0-e2*SinPhi*SinPhi);
103 dZ_new = N*e2*SinPhi;
104 if ( Math.abs(dZ-dZ_new) < epsRequ ) break;
105 dZ = dZ_new;
106 }
107
108 // Longitude, latitude, altitude
109 //double[] LLA = new double[3];
110 LLA[1] = Math.atan2( Y, X ); // longitude, lon
111 LLA[0] = Math.atan2( ZdZ, Math.sqrt(rho2) ); // latitude, lat
112 LLA[2] = Nh - N; // altitute, h
113
114 //System.out.println("LLA[1]: "+ LLA[1]);
115 //LLA[1] = LLA[1] -(280.4606 +360.9856473*d)*Math.PI/180.0; // shift based on time
116 // add fidelity to the line above
117 LLA[1] = LLA[1] - earthRotationDeg(d)*Math.PI/180.0; // shift based on time
118 double div = Math.floor(LLA[1]/(2*Math.PI));
119 LLA[1] = LLA[1] - div*2*Math.PI;
120 if(LLA[1] > Math.PI)
121 {
122 LLA[1] = LLA[1]- 2.0*Math.PI;
123 }
124
125 //System.out.println("LLA[1]a: "+ LLA[1]);
126
127 return LLA; //h
128
129 } // calculateGeodeticLLA
130
131 // SEG 10 June 2009 - help standardize earth Rotations
132 private static double earthRotationDeg(double d) // days since y2K
133 {
134 // LLA[1] = LLA[1] -(280.4606 +360.9856473*d)*Math.PI/180.0; // shift based on time
135
136 // calculate T
137 double T = (d)/36525.0;
138
139 // do calculation
140 return ( (280.46061837 + 360.98564736629*(d)) + 0.000387933*T*T - T*T*T/38710000.0) % 360.0;
141
142 } // earthRotationRad
143
144 /**
145 * calculate the pointing information Azumuth, Elevation, and Range (AER) to
146 * a satellite from a location on Earth (given Lat, Long, Alt)
147 * if elevation >=0 then sat is above local horizon
148 * @param currentJulianDate Julian Date for AER calculation (corresponds to ECI position)
149 * @param lla_deg_m lat long and alt of station in deg/deg/meters (Geodetic)
150 * @param eci_pos ECI position of object in meters (sat)
151 * @return Azumuth [deg], Elevation [deg], and Range vector [m]
152 */
153 public static double[] calculate_AER(double currentJulianDate,double[] lla_deg_m, double[] eci_pos)
154 { double[] aer = new double[3];
155
156 // 0th step get local mean Sidereal time
157 // first get mean sidereal time for this station - since we use it twice
158 double thetaDeg = Sidereal.Mean_Sidereal_Deg(currentJulianDate-AstroConst.JDminusMJD, lla_deg_m[1]);
159 // first calculate ECI position of Station
160 double[] eciGS = calculateECIposition(lla_deg_m,thetaDeg);
161
162 // find the vector between pos and GS
163 double[] rECI = MathUtils.sub(eci_pos, eciGS);
164
165 // calculate range
166 aer[2] = MathUtils.norm(rECI);
167
168 // now transform ECI to topocentric-horizon system (SEZ) (use Geodetic Lat, not geocentric)
169 double[] rSEZ = eci2sez(rECI,thetaDeg,lla_deg_m[0]); // ECI vec, sidereal in Deg, latitude in deg
170
171 // compute azimuth [radians] -> Deg
172 //aer[0] = Math.atan(-rSEZ[1]/rSEZ[0]) * 180.0/Math.PI;
173 aer[0] = Math.atan2(-rSEZ[0], rSEZ[1]) * 180.0/Math.PI;
174
175 //System.out.println("aer[0]_0=" + aer[0] + ", rSEZ[-0,1]=" + (-rSEZ[0]) + ", " +rSEZ[1]);
176
177 // do conversions so N=0, S=180, NW=270
178 if(aer[0] <= 0)
179 {
180 aer[0] = Math.abs(aer[0]) + 90;
181 }
182 else
183 {
184 if(aer[0]<= 90) //(between 0 and 90)
185 {
186 aer[0] = -1.0*aer[0] + 90.0;
187 }
188 else // between 90 and 180
189 {
190 aer[0] = -1.0*aer[0] + 450.0;
191 }
192 }
193
194 // compute elevation [radians]
195 aer[1] = Math.asin(rSEZ[2] / aer[2]) * 180.0/Math.PI;
196
197 //System.out.println("SEZ: " + rSEZ[0] + ", " + rSEZ[1] + ", " + rSEZ[2]);
198
199 return aer;
200 } // calculate_AER
201
202 /**
203 * Calculate ECI position from local mean sidereal time and geodetic lat long alt
204 * @param lla_deg_m lat long and alt of station in deg/deg/meters (Geodetic)
205 * @param theta local mean sidereal time (Degrees)
206 * @return ECI position (meters)
207 */
208 public static double[] calculateECIposition(double[] lla_deg_m, double theta)
209 {
210 // calculate the ECI j2k position vector of the ground station at the current time
211 double [] eciVec = new double[3];
212
213 // // calculate geocentric latitude - using non spherical earth (in radians)
214 // // http://celestrak.com/columns/v02n03/
215 // double geocentricLat = Math.atan( Math.pow(1.0-AstroConst.f_Earth, 2.0) * Math.tan( lla_deg_m[0]*Math.PI/180.0 ) ); // (1-f)^2 tan(?).
216 //
217 // eciVec[2] = AstroConst.R_Earth * Math.sin( geocentricLat ); //lla_deg_m[0]*Math.PI/180.0 );
218 // double r = AstroConst.R_Earth * Math.cos( geocentricLat ); //lla_deg_m[0]*Math.PI/180.0 );
219 // eciVec[0] = r * Math.cos(theta*Math.PI/180.0);
220 // eciVec[1] = r * Math.sin(theta*Math.PI/180.0);
221
222 // alternate way to calcuate ECI position - using earth flattening
223 // http://celestrak.com/columns/v02n03/
224 double C = 1.0 / Math.sqrt( 1.0+AstroConst.f_Earth*(AstroConst.f_Earth-2.0)*Math.pow(Math.sin(lla_deg_m[0]*Math.PI/180.0 ),2.0) );
225 double S = Math.pow(1.0-AstroConst.f_Earth, 2.0) * C;
226
227 eciVec[0] = AstroConst.R_Earth * C * Math.cos(lla_deg_m[0]*Math.PI/180.0)*Math.cos(theta*Math.PI/180.0);
228 eciVec[1] = AstroConst.R_Earth * C * Math.cos(lla_deg_m[0]*Math.PI/180.0)*Math.sin(theta*Math.PI/180.0);
229 eciVec[2] = AstroConst.R_Earth * S * Math.sin(lla_deg_m[0]*Math.PI/180.0);
230
231 return eciVec;
232
233 } //calculateECIposition
234
235 /**
236 * transform ECI to topocentric-horizon system (SEZ) (south-East-Zenith)
237 * @param rECI position in ECI coordinates (meters)
238 * @param thetaDeg local sidereal time (degrees)
239 * @param latDeg observer's latitude (degrees)
240 * @return topocentric-horizon system (SEZ) (south-East-Zenith)
241 */
242 public static double[] eci2sez(double[] rECI,double thetaDeg,double latDeg)
243 {
244 double[] rSEZ = new double[3]; // new postion in SEZ coorinates
245
246 //? (the local sidereal time) -> (thetaDeg*Math.PI)
247 //? (the observer's latitude) - > (latDeg*Math.PI)
248 rSEZ[0] = Math.sin(latDeg*Math.PI/180.0) * Math.cos(thetaDeg*Math.PI/180.0) * rECI[0] + Math.sin(latDeg*Math.PI/180.0) * Math.sin(thetaDeg*Math.PI/180.0) * rECI[1] - Math.cos(latDeg*Math.PI/180.0) * rECI[2];
249 rSEZ[1] = -Math.sin(thetaDeg*Math.PI/180.0) * rECI[0] + Math.cos(thetaDeg*Math.PI/180.0) * rECI[1]; rSEZ[2] = Math.cos(latDeg*Math.PI/180.0) * Math.cos(thetaDeg*Math.PI/180.0) * rECI[0] +
250 Math.cos(latDeg*Math.PI/180.0) * Math.sin(thetaDeg*Math.PI/180.0) * rECI[1] + Math.sin(latDeg*Math.PI/180.0) * rECI[2];
251
252 return rSEZ;
253 }
254 }