| 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; |
