| 1 |
greg |
1.1 |
{
|
| 2 |
|
|
Coordinate mapping for point above Earth
|
| 3 |
|
|
assuming uniform projection of latitude and longitude
|
| 4 |
|
|
|
| 5 |
|
|
5/1/96 Greg Ward
|
| 6 |
|
|
|
| 7 |
|
|
A1 = latitude position (degrees)
|
| 8 |
|
|
A2 = longitude position (degrees)
|
| 9 |
|
|
A3 = height above Earth (kilometers)
|
| 10 |
|
|
A4 = Julian date (1-365)
|
| 11 |
|
|
A5 = Greenwich Mean Time (0-24)
|
| 12 |
|
|
|
| 13 |
|
|
By modifying this file so that the above arguments are defined
|
| 14 |
|
|
instead as constants, significant gains in speed are possible.
|
| 15 |
|
|
(Be sure to change Su, Sv and Sw to constant definitions also.)
|
| 16 |
|
|
|
| 17 |
|
|
Y-axis points North, Z-axis points up (away from Earth)
|
| 18 |
|
|
}
|
| 19 |
|
|
{
|
| 20 |
|
|
Example usage:
|
| 21 |
|
|
|
| 22 |
|
|
void colorpict epict
|
| 23 |
|
|
7 Eval Eval Eval Earth.pic Earth.cal u v
|
| 24 |
|
|
0
|
| 25 |
|
|
5 lat long height JDate GMT
|
| 26 |
|
|
|
| 27 |
|
|
epict glow eglow
|
| 28 |
|
|
0
|
| 29 |
|
|
0
|
| 30 |
|
|
4 1 1 1 0
|
| 31 |
|
|
|
| 32 |
|
|
eglow source earth
|
| 33 |
|
|
0
|
| 34 |
|
|
0
|
| 35 |
|
|
4 0 0 -1 max_angle
|
| 36 |
|
|
|
| 37 |
|
|
where:
|
| 38 |
|
|
max_angle = 2*180/PI * asin(EarthRad/(EarthRad+height))
|
| 39 |
|
|
|
| 40 |
|
|
If you are really close to the Earth (below atmosphere), use:
|
| 41 |
|
|
|
| 42 |
|
|
!gensky month day time -a lat -o lon -m std_meridian
|
| 43 |
|
|
|
| 44 |
|
|
skybright colorpict epict
|
| 45 |
|
|
7 noop noop noop Earth.pic Earth.cal u v
|
| 46 |
|
|
0
|
| 47 |
|
|
3 lat long height
|
| 48 |
|
|
|
| 49 |
|
|
epict glow eglow
|
| 50 |
|
|
0
|
| 51 |
|
|
0
|
| 52 |
|
|
4 1 1 1 0
|
| 53 |
|
|
|
| 54 |
|
|
eglow source earth
|
| 55 |
|
|
0
|
| 56 |
|
|
0
|
| 57 |
|
|
4 0 0 -1 max_angle
|
| 58 |
|
|
}
|
| 59 |
|
|
|
| 60 |
|
|
EarthRad : 6378.5; { Earth radius in kilometers }
|
| 61 |
|
|
EarthTilt : 23.45*DEGREE; { Earth's orbital axis tilt }
|
| 62 |
|
|
EarthIrrad : 770 {watts/meter^2}; { Maximum Earth irradiance }
|
| 63 |
|
|
PA : 256/512; { Picture Height/Width }
|
| 64 |
|
|
|
| 65 |
|
|
{ Change the following to constants to speed things up }
|
| 66 |
|
|
Lat = arg(1)*DEGREE; { Position latitude }
|
| 67 |
|
|
Lon = arg(2)*DEGREE; { Position longitude }
|
| 68 |
|
|
Height = arg(3); { Position height above Earth (km) }
|
| 69 |
|
|
JDate = arg(4); { Julian date (1-365) }
|
| 70 |
|
|
GMT = arg(5); { Greenwich Mean Time (0-24) }
|
| 71 |
|
|
|
| 72 |
|
|
GOu = EarthRad*sin(Lon)*cos(Lat); { Starting Earth position }
|
| 73 |
|
|
GOv = EarthRad*cos(Lon)*cos(Lat);
|
| 74 |
|
|
GOw = EarthRad*sin(Lat);
|
| 75 |
|
|
|
| 76 |
|
|
raduv = max(FTINY, sqrt(GOu*GOu + GOv*GOv)); { Earthwise direction vector }
|
| 77 |
|
|
GVXu = -GOv / raduv;
|
| 78 |
|
|
GVXv = GOu / raduv;
|
| 79 |
|
|
GVXw = 0;
|
| 80 |
|
|
GVYu = (-GOw*GOu) / EarthRad / raduv;
|
| 81 |
|
|
GVYv = (GOw*-GOv) / EarthRad / raduv;
|
| 82 |
|
|
GVYw = (GOu*GOu + GOv*GOv) / EarthRad / raduv;
|
| 83 |
|
|
radxy = max(FTINY, sqrt(Dx*Dx + Dy*Dy));
|
| 84 |
|
|
GDu = (GVXu*Dx + GVYu*Dy) / radxy;
|
| 85 |
|
|
GDv = (GVXv*Dx + GVYv*Dy) / radxy;
|
| 86 |
|
|
GDw = (GVXw*Dx + GVYw*Dy) / radxy;
|
| 87 |
|
|
|
| 88 |
|
|
tana = if(FTINY+Dz, radxy/FTINY, radxy/-Dz);
|
| 89 |
|
|
|
| 90 |
|
|
det = sq(EarthRad+Height) - (tana*tana + 1)*(2*EarthRad*Height + Height*Height);
|
| 91 |
|
|
|
| 92 |
|
|
a = if(tana-FTINY, b/tana, Height);
|
| 93 |
|
|
|
| 94 |
|
|
b = if(tana-FTINY, (EarthRad + Height - Sqrt(det)) / (tana + 1/tana), 0);
|
| 95 |
|
|
|
| 96 |
|
|
c = EarthRad + Height - a;
|
| 97 |
|
|
|
| 98 |
|
|
digfactor = 1 - (a - Height)/EarthRad;
|
| 99 |
|
|
{ Global position vector }
|
| 100 |
|
|
GPu = GOu*digfactor + GDu*b;
|
| 101 |
|
|
GPv = GOv*digfactor + GDv*b;
|
| 102 |
|
|
GPw = GOw*digfactor + GDw*b;
|
| 103 |
|
|
{ Earthly coordinates }
|
| 104 |
|
|
longitude = atan2(GPu, GPv);
|
| 105 |
|
|
latitude = Asin(GPw/EarthRad);
|
| 106 |
|
|
{ Image lookup coordinates }
|
| 107 |
|
|
u = if(PA-1, 1, 1/PA) * (.5 - longitude/(2*PI));
|
| 108 |
|
|
v = if(PA-1, PA, 1) * (.5 + latitude/PI);
|
| 109 |
|
|
{ Solar declination (lat.) }
|
| 110 |
|
|
Slat = EarthTilt * sin(2*PI/368 * (JDate - 81));
|
| 111 |
|
|
{ Solar longitude }
|
| 112 |
|
|
Slon = 2*PI/24 * (GMT - 12);
|
| 113 |
|
|
{ Solar position vector }
|
| 114 |
|
|
Su = EarthRad*sin(Slon)*cos(Slat);
|
| 115 |
|
|
Sv = EarthRad*cos(Slon)*cos(Slat);
|
| 116 |
|
|
Sw = EarthRad*sin(Slat);
|
| 117 |
|
|
{ Solar cosine }
|
| 118 |
|
|
Scos = (Su*GPu + Sv*GPv + Sw*GPw)/(EarthRad*EarthRad);
|
| 119 |
|
|
{ Earth radiance multiplier }
|
| 120 |
|
|
Eradmult = EarthIrrad / PI * noneg(Scos);
|
| 121 |
|
|
{ Earth radiance }
|
| 122 |
|
|
Eval(rho) = Eradmult * rho;
|
