001 package edu.nrao.sss.astronomy;
002
003 import java.math.BigDecimal;
004 import java.util.Date;
005
006 import static edu.nrao.sss.astronomy.CelestialCoordinateSystem.*;
007
008 import edu.nrao.sss.geom.EarthPosition;
009 import edu.nrao.sss.measure.Angle;
010 import edu.nrao.sss.measure.ArcUnits;
011 import edu.nrao.sss.measure.Latitude;
012 import edu.nrao.sss.measure.LocalSiderealTime;
013 import edu.nrao.sss.measure.Longitude;
014
015 /**
016 * A converter that converts only equatorial positions to or from
017 * horizontal positions.
018 * <p>
019 * <b>Version Info:</b>
020 * <table style="margin-left:2em">
021 * <tr><td>$Revision: 1242 $</td></tr>
022 * <tr><td>$Date: 2008-04-28 08:46:17 -0600 (Mon, 28 Apr 2008) $</td></tr>
023 * <tr><td>$Author: dharland $ (last person to modify)</td></tr>
024 * </table></p>
025 *
026 * @author David M. Harland
027 * @since 2008-02-27
028 */
029 public class EquatorialHorizontalConverter
030 implements CelestialCoordinateConverter
031 {
032 public SkyPosition createFrom(SkyPosition position,
033 CelestialCoordinateSystem toSystem,
034 Epoch toEpoch,
035 EarthPosition observer,
036 LocalSiderealTime lst)
037 throws CoordinateConversionException
038 {
039 //TODO Do we implicitly support only J2000?
040 // If so, throw exception for B1950?
041
042 CelestialCoordinateSystem fromSystem = position.getCoordinateSystem();
043
044 SkyPosition sp = null;
045
046 //TODO do null-conversion as a courtesy?
047 // Note: non-trivial if HORZ to HORZ and position has
048 // changed.
049
050 if (fromSystem.equals(EQUATORIAL) && toSystem.equals(HORIZONTAL))
051 {
052 sp = equatorialToHorizontal(position, observer, lst);
053 }
054 else if (fromSystem.equals(HORIZONTAL) && toSystem.equals(EQUATORIAL))
055 {
056 sp = horizontalToEquatorial(position, observer, lst);
057 }
058 else
059 {
060 //TODO throw excep
061 }
062
063 return sp;
064 }
065
066 /** Creates an equatorial position from a horizontal. */
067 private SkyPosition
068 horizontalToEquatorial(SkyPosition horizPosition,
069 EarthPosition observerPosition,
070 LocalSiderealTime lst)
071 {
072 Date time = lst.toDate();
073
074 Angle elAngle = horizPosition.getLatitude(time).toAngle();
075 Angle azAngle = horizPosition.getLongitude(time).toAngle();
076 Angle latAngle = observerPosition.getLatitude().toAngle();
077
078 //This algorithm is one that avoids arcsin and arccosine.
079 //It comes from http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?dh2e.src
080 //(SLALIB) and works by transforming rectangular to polar coordinates.
081
082 BigDecimal SA = new BigDecimal(azAngle.sine());
083 BigDecimal CA = new BigDecimal(azAngle.cosine());
084 BigDecimal SE = new BigDecimal(elAngle.sine());
085 BigDecimal CE = new BigDecimal(elAngle.cosine());
086 BigDecimal SP = new BigDecimal(latAngle.sine());
087 BigDecimal CP = new BigDecimal(latAngle.cosine());
088
089 BigDecimal X = CA.negate().multiply(CE).multiply(SP);
090 BigDecimal temp = SE.multiply(CP);
091 X = X.add(temp);
092
093 BigDecimal Y = SA.negate().multiply(CE);
094
095 BigDecimal Z = CA.multiply(CE).multiply(CP);
096 temp = SE.multiply(SP);
097 Z = Z.add(temp);
098
099 BigDecimal Rsq = X.multiply(X);
100 temp = Y.multiply(Y);
101 Rsq = Rsq.add(temp);
102
103 double R = Math.sqrt(Rsq.doubleValue());
104
105 double HA = (R == 0.0) ? 0.0: Math.atan2(Y.doubleValue(), X.doubleValue());
106
107 Angle hourAngle = new Angle(BigDecimal.valueOf(HA), ArcUnits.RADIAN);
108
109 Angle declination = new Angle(BigDecimal.valueOf(Math.atan2(Z.doubleValue(), R)),
110 ArcUnits.RADIAN);
111
112 Longitude rightAscension = lst.getRightAscension(new Longitude(hourAngle));
113
114 //Construct the new position
115 SimpleSkyPosition eqPos =
116 new SimpleSkyPosition(CelestialCoordinateSystem.EQUATORIAL,
117 Epoch.J2000);
118 eqPos.setLongitude(rightAscension);
119 eqPos.setLatitude(new Latitude(declination));
120
121 return eqPos;
122 }
123
124 /** Creates a horizontal position from an equatorial. */
125 private SkyPosition
126 equatorialToHorizontal(SkyPosition eqPosition,
127 EarthPosition observerPosition,
128 LocalSiderealTime lst)
129 {
130 Date time = lst.toDate();
131
132 Angle decAngle = eqPosition.getLatitude(time).toAngle();
133 Angle latAngle = observerPosition.getLatitude().toAngle();
134 Angle hourAngle = lst.getHourAngle(eqPosition.getLongitude(time)).toAngle();
135
136 //This algorithm is one that avoids arcsin and arccosine.
137 //It comes from http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?de2h.src
138 //(SLALIB) and works by transforming rectangular to polar coordinates.
139
140 BigDecimal SH = new BigDecimal(hourAngle.sine());
141 BigDecimal CH = new BigDecimal(hourAngle.cosine());
142 BigDecimal SD = new BigDecimal(decAngle.sine());
143 BigDecimal CD = new BigDecimal(decAngle.cosine());
144 BigDecimal SP = new BigDecimal(latAngle.sine());
145 BigDecimal CP = new BigDecimal(latAngle.cosine());
146
147 BigDecimal X = CH.negate().multiply(CD).multiply(SP);
148 BigDecimal temp = SD.multiply(CP);
149 X = X.add(temp);
150
151 BigDecimal Y = SH.negate().multiply(CD);
152
153 BigDecimal Z = CH.multiply(CD).multiply(CP);
154 temp = SD.multiply(SP);
155 Z = Z.add(temp);
156
157 BigDecimal Rsq = X.multiply(X);
158 temp = Y.multiply(Y);
159 Rsq = Rsq.add(temp);
160
161 double R = Math.sqrt(Rsq.doubleValue());
162
163 double A = (R == 0.0) ? 0.0: Math.atan2(Y.doubleValue(), X.doubleValue());
164
165 Angle azimuth = new Angle(BigDecimal.valueOf(A), ArcUnits.RADIAN);
166
167 Angle elevation = new Angle(BigDecimal.valueOf(Math.atan2(Z.doubleValue(), R)),
168 ArcUnits.RADIAN);
169
170 //Construct the new position
171 SimpleSkyPosition horizPos =
172 new SimpleSkyPosition(CelestialCoordinateSystem.HORIZONTAL,
173 eqPosition.getEpoch());
174 horizPos.setLongitude(new Longitude(azimuth));
175 horizPos.setLatitude(new Latitude(elevation));
176
177 return horizPos;
178 }
179
180 //============================================================================
181 //
182 //============================================================================
183 /*
184 public static void main(String... args)
185 {
186 LocalSiderealTime lst = new LocalSiderealTime();
187 System.out.println("LST = "+lst.toDate()+"\n");
188
189 EarthPosition observer = new EarthPosition(lst.getLocation(), Latitude.parse("34d 04' 43.497\""), lst.getTimeZone());
190 System.out.println("Longitude = "+observer.getLongitude().toStringDms()+
191 ", Latitude = "+observer.getLatitude().toStringDms()+"\n");
192
193 SimpleSkyPosition eqPos = new SimpleSkyPosition(CelestialCoordinateSystem.EQUATORIAL, Epoch.J2000);
194 System.out.print("Orig RA/Dec: " ); displayPosition(eqPos);
195
196 EquatorialHorizontalConverter converter = new EquatorialHorizontalConverter();
197 SkyPosition horzPos = converter.equatorialToHorizontal(eqPos, observer, lst);
198 System.out.print("Conv Az/El: " ); displayPosition(horzPos);
199
200 SkyPosition eqPos2 = converter.horizontalToEquatorial(horzPos, observer, lst);
201 System.out.print("Converted back: " ); displayPosition(eqPos2);
202
203 //The elevation of the celestial pole should equal latitude
204 System.out.println();
205 eqPos.setLatitude(new Latitude(90.0));
206 eqPos.setLongitude(new Longitude(0.0));
207 horzPos = converter.equatorialToHorizontal(eqPos, observer, lst);
208 displayPosition(horzPos);
209
210 eqPos2 = converter.horizontalToEquatorial(horzPos, observer, lst);
211 displayPosition(eqPos2);
212
213 System.out.println();
214 eqPos.setLongitude(new Longitude(60.0));
215 eqPos.setLatitude(new Latitude(50.0));
216 System.out.print("Orig RA/Dec: " ); displayPosition(eqPos);
217 horzPos = converter.equatorialToHorizontal(eqPos, observer, lst);
218 System.out.print("Conv Az/El: " ); displayPosition(horzPos);
219 eqPos2 = converter.horizontalToEquatorial(horzPos, observer, lst);
220 System.out.print("Converted back: " ); displayPosition(eqPos2);
221 }
222 private static void displayPosition(SkyPosition pos)
223 {
224 System.out.print(pos.getCoordinateSystem());
225 if (pos.getCoordinateSystem().equals(EQUATORIAL))
226 System.out.print("-"+pos.getEpoch());
227 System.out.print(", ");
228 System.out.print(pos.getLongitude().toUnits(ArcUnits.DEGREE)); //.toStringHms());
229 System.out.print(", ");
230 System.out.print(pos.getLatitude().toUnits(ArcUnits.DEGREE)); //.toStringDms());
231 System.out.println();
232 }
233 */
234 }