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}