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