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