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