Attachment 'll288.c'
Download 1 // ll288.c
2 // compute McIlwain L layer versus angle in server sky orbit
3 // also compute the histogram of probability versus L shell
4 //
5 // gcc -o ll288 ll288.c -lm
6
7 #define NAME "M288"
8
9 // #define NPOINTS 37 //
10 // #define NBINS 10 //
11 // #define NSAMPLE 20 //
12
13 #define NPOINTS 361 //
14 #define NBINS 100 //
15 #define NSAMPLE 360000 //
16
17 #define DEBUG 1 //
18 #define AMIN -180.0 // orbit longitude
19 #define AMAX 180.0 //
20
21 #define RE 6371.0 // McIlwain earth radius kilometers
22 #define M288 12789.0 // M288
23 #define M360 14441.0 // M360
24 #define M480 16756.0 // M480
25 #define RORB (M288) // Orbit radius kilometers
26
27 #define ANG_N 10.26 // field tilt from north
28 #define ANG_E -71.78 // field tilt to east
29
30 #define OFF_R 552.0 // field offset in kilometers
31 #define OFF_N 22.2 // field offset tilt from north
32 #define OFF_E 141.6 // field offset tilt towards east
33
34 #define RAD (3.1415926536/180.0)
35
36 #include <math.h>
37 #include <string.h>
38 #include <stdio.h>
39 #include <stdlib.h>
40
41 //-----------------------------------------------------------------
42 double h_o ;
43 double x_o ; // equatorial plane
44 double y_o ; // NORTH (as POVRAY)
45 double z_o ; // equatorial plane (as POVRAY)
46
47 double lm, bm, phim ;
48
49 double lmin, lmax ;
50
51 //-----------------------------------------------------------------
52 void lshell( double ang ) {
53
54 // position in earth radius=1 coordinates
55
56 double rorb = RORB/RE ;
57 double xorb = -rorb * sin( RAD * ang );
58 double yorb = 0.0 ;
59 double zorb = -rorb * cos( RAD * ang );
60
61 // position relative to dipole offset
62
63 double x1 = xorb - x_o ;
64 double y1 = yorb - y_o ;
65 double z1 = zorb - z_o ;
66
67 // rotate relative to dipole east-west
68
69 double a2 = RAD * ANG_E ;
70 double x2 = x1 * cos(a2) + z1 * sin(a2) ;
71 double y2 = y1 ;
72 double z2 = z1 * cos(a2) - x1 * sin(a2) ;
73
74 // rotate relative to dipole north-south
75
76 double am = RAD * ANG_N ;
77 double xm = x2 * cos(am) + y2 * sin(am) ;
78 double ym = y2 * cos(am) - x2 * sin(am) ;
79 double zm = z2 ;
80 double xzm = sqrt( xm*xm + zm*zm ) ;
81 double rm = sqrt( xzm*xzm + ym*ym ) ;
82
83 // find McIlwain parameters
84
85 phim = atan2( ym , xzm ) ;
86 double sinm = xzm / rm ;
87 double cosm = ym / rm ;
88 lm = rm / ( sinm * sinm ) ;
89 bm = sqrt( 3.0*cosm*cosm+1.0 )
90 / ( rm*rm*rm ) ;
91
92 /*
93 printf( "%10.3f_ang%10.3f_xm%10.3f_ym%10.3f_zm%10.3f_phim\n",
94 ang, xm, ym, zm, phim );
95 */
96 }
97
98 //-----------------------------------------------------------------
99 int main() {
100
101 FILE *datfile ; //
102 char gnuplot[200] ; //
103 char filename[80] ; //
104 int iang ; //
105
106 lmin = 1000.0 ;
107 lmax = -1000.0 ;
108
109 h_o = OFF_R * cos( RAD * OFF_N ) ;
110 x_o = h_o * sin( RAD * OFF_E ) / RE ;
111 y_o = OFF_R * sin( RAD * OFF_N ) / RE ;
112 z_o = -h_o * cos( RAD * OFF_E ) / RE ;
113
114 // --------------------------------------------------------------------
115 // compute angle loop
116
117 sprintf( filename, "%ss.dat", NAME );
118 datfile = fopen( filename, "wb" );
119 for( iang=0 ; iang < NPOINTS ; iang++ ) {
120 double ang = AMIN + (AMAX-AMIN)*((double)iang )/((double)(NPOINTS-1));
121
122 lshell( ang ) ;
123 if( lmin > lm ) { lmin = lm ; }
124 if( lmax < lm ) { lmax = lm ; }
125
126 double lamm = fabs( phim/RAD ) ;
127
128 fprintf( datfile, "%16.4f%16.4f%16.4f%16.4f\n", ang, lm, bm, lamm ) ;
129 }
130 fclose( datfile );
131
132 lmax += 1e-4 ;
133 lmin -= 1e-4 ;
134
135 // histogram
136
137 double prob[NBINS], lval ; //
138 double delta = ((double) NBINS )/((double) NSAMPLE );
139
140 // histogram initialization
141 for( iang=0 ; iang < NBINS ; iang++ ) {
142 prob[ iang ] = 0.0 ;
143 }
144
145 double lscale = ((double)NBINS) /(lmax-lmin);
146
147 // histogram compute
148 for( iang=0 ; iang < NSAMPLE ; iang++ ) {
149 double ang = AMIN + (AMAX-AMIN)*((double)iang )/((double)(NSAMPLE));
150 lshell( ang ) ;
151 int bin = (int) ( (lm-lmin)*lscale ) ;
152 prob[ bin ] += delta ;
153 }
154
155 // histogram printout
156 sprintf( filename, "%sh.dat", NAME );
157 datfile = fopen( filename, "wb" );
158
159 double dlval = (lmax-lmin)/((double) NBINS );
160
161 lval = lmin + 0.5*dlval ;
162 fprintf( datfile, "%16.4f%16.4f%16.4f\n", 0.0, lval ) ;
163
164 for( iang=0 ; iang < NBINS ; iang++ ) {
165 lval = lmin + dlval * (0.5+(double)iang) ;
166 fprintf( datfile, "%16.4f%16.4f%16.4f\n", prob[ iang ], lval ) ;
167 }
168
169 lval = lmax - 0.5*dlval ;
170 fprintf( datfile, "%16.4f%16.4f%16.4f\n", 0.0, lval ) ;
171
172 fclose( datfile );
173
174 // prepare gnuplot controp file
175 sprintf( gnuplot, "cat XXXXX.gp | sed -e \"s/XXXXX/%s/\" > %s.gp",
176 NAME, NAME );
177 system( gnuplot );
178
179 // run gnuplot
180 sprintf( gnuplot, "/usr/local/bin/gp %s.gp", NAME );
181 system( gnuplot );
182 return(0);
183 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.