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