001/* 002 * $Id: GeoFunctions.java,v 1.4 2011/03/24 16:06:33 davep Exp $ 003 * 004 * This file is part of McIDAS-V 005 * 006 * Copyright 2007-2011 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 */ 030package edu.wisc.ssec.mcidasv.data.adde.sgp4; 031 032/** 033 * Earth computation functions 034 */ 035public 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] + 252Math.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}