00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070 #include "wsa_pkg.h"
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080 float wsa_f(float f, float th)
00081 {
00082 float fs_exp, bb, kk, th0, th_exp, all_exp;
00083 float fac0, expo, fac1, spd;
00084 fs_exp = 0.22; bb = 6.08; kk = 3.30;
00085 th0 = 2.50; th_exp = 1.50; all_exp = 3.40;
00086 fac0 = pow((1. + fabs(f)), fs_exp);
00087 expo = 1 - pow((th / th0), th_exp);
00088 fac1 = bb - kk * exp(expo) / exp(1.0);
00089 spd = 265. + 1.5 * pow(fac1, all_exp) / fac0;
00090 return spd;
00091 }
00092
00093
00094
00095
00096
00097
00098
00099 float gcd(float lon1, float lat1, float lon2, float lat2)
00100 {
00101 double dlon = lon1 - lon2;
00102 double dlat = lat1 - lat2;
00103 if ((fabs(dlon) < 1.E-5) && (fabs(dlat) < 1.E-5))
00104 return 0.;
00105 double num1 = cos(lat2) * sin(dlon);
00106 double num2 = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dlon);
00107 double num = sqrt(num1 * num1 + num2 * num2);
00108 double den = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(dlon);
00109 float result = fabs(atan2(num, den));
00110 return (result);
00111 }
00112
00113
00114
00115 void smoothfp(int *fpmap, int np, int nt)
00116 {
00117 int Dpt, Dpt0;
00118 for (int j = 0; j < nt; j++) {
00119 for (int i = 0; i < np; i++) {
00120 Dpt = j * np + i;
00121 if (fpmap[Dpt]) continue;
00122 int count = 0;
00123 for (int ii = -1; ii <= 1; ii++) {
00124 for (int jj = -1; jj <= 1; jj++) {
00125 if (ii == 0 && jj == 0) continue;
00126 int lon = i + ii;
00127 int lat = j + jj;
00128 if (lon < 0) lon += np;
00129 if (lon >= np) lon -= np;
00130 if (lat < 0) {
00131 lat = abs(lat);
00132 lon = lon + np / 2;
00133 if (lon >= np) lon -= np;
00134 }
00135 if (lat >= nt) {
00136 lat = 2 * nt - lat - 1;
00137 lon = lon + np / 2;
00138 if (lon >= np) lon -= np;
00139 }
00140 Dpt0 = lat * np + lon;
00141 if (fpmap[Dpt0]) count++;
00142 }
00143 }
00144 if (count >= 5) fpmap[Dpt] = 1;
00145 }
00146 }
00147 }
00148
00149
00150
00151
00152
00153 void getbd(int *fpmap, float *fte, float *fpph, float *fpth,
00154 int np, int nt)
00155 {
00156 float dph = 360. / np, dth = 180. / nt;
00157
00158 int Dpt, Dpt0;
00159 for (int j = 0; j < nt; j++) {
00160 for (int i = 0; i < np; i++) {
00161 Dpt = j * np + i;
00162 if (fabs(fte[Dpt]) > 1.E-5) {
00163 int ii = (int)(fpph[Dpt] / DTOR / dph);
00164 int jj = nt - 1 - (int)(fpth[Dpt] / DTOR / dth);
00165 Dpt0 = jj * np + ii;
00166 fpmap[Dpt0] = 1;
00167 }
00168 }
00169 }
00170
00171 for (int i = 0; i < 15; i++)
00172 smoothfp(fpmap, np, nt);
00173 }
00174
00175
00176
00177
00178
00179
00180 float ang_sep(float ph, float th, int *fpmap, int np, int nt)
00181 {
00182 float dph = 360. / np, dth = 180. / nt;
00183 float xx = ph, yy = M_PI / 2. - th;
00184 float dist = 181.0, dist_t = 0.0;
00185
00186 int i, ii, jj;
00187 int Dpt;
00188 float dx, dy;
00189 int x_grid, y_grid;
00190 int box_x, box_y;
00191 int left, right, up, down;
00192 float lon, lat;
00193
00194 dx = dph * DTOR;
00195 dy = dth * DTOR;
00196 x_grid = (int)(ph / dx);
00197 y_grid = (int)((M_PI - th) / dy);
00198 box_x = (int)(BOX / dph / cos(yy));
00199 box_y = (int)(BOX / dth);
00200 if (box_x >= (np / 2)) {
00201 left = 0;
00202 right = np - 1;
00203 } else {
00204 left = x_grid - box_x;
00205 right = x_grid + box_x;
00206 }
00207 up = y_grid + box_y;
00208 down = y_grid - box_y;
00209 for (int aa = left; aa <= right; aa++) {
00210 if (aa < 0)
00211 i = aa + np;
00212 else if (aa >= np)
00213 i = aa - np;
00214 else
00215 i = aa;
00216 for (int bb = down; bb <= up; bb++) {
00217 if (bb < 0) {
00218 ii = i + np / 2;
00219 if (ii >= np) ii -= np;
00220 jj = abs(bb);
00221 } else if (bb >= nt) {
00222 ii = i + np / 2;
00223 if (ii >= np) ii -= np;
00224 jj = 2 * nt - bb - 1;
00225 } else {
00226 ii = i;
00227 jj = bb;
00228 }
00229 Dpt = jj * np + ii;
00230 if ((ii != x_grid) && (jj != y_grid) && !fpmap[Dpt]) {
00231 lon = (ii + 0.5) / np * 2 * M_PI;
00232 lat = (jj + 0.5) / nt * M_PI - M_PI / 2.;
00233 dist_t = gcd(xx, yy, lon, lat) / DTOR;
00234 if (dist_t < dist) dist = dist_t;
00235 }
00236 }
00237 }
00238 return dist;
00239 }
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 void fte_global(float *g, float *h, struct Grid *grid, float *fte,
00250 float *fpph, float *fpth, float *brss, int lmax,
00251 void (*integrator)(struct Point *, struct Point *, float *, float *, int, float, float,
00252 void (*)(float *, float *, struct Point *, float *, int, float)))
00253 {
00254 int np = grid->np, nt = grid->nt;
00255 int failed = 0;
00256
00257 struct Point *pt0 = (struct Point *)(malloc(sizeof(struct Point)));
00258 pt0->r = RS;
00259 struct FldLn *fl = (struct FldLn *)(malloc(sizeof(struct FldLn)));
00260 fl->num = 0;
00261 void (*Bfunc)(float *, float *, struct Point *, float *, int, float);
00262 Bfunc = pfss0;
00263
00264 for (int i = 0; i < np; i++) {
00265 for (int j = 0; j < nt; j++) {
00266 pt0->t = grid->th[j];
00267 pt0->p = grid->ph[i];
00268
00269 trace(g, h, pt0, fl, lmax, integrator, Bfunc, RS, R0 + 0.001, 0.0);
00270
00271 int Dpt = j * np + i;
00272 brss[Dpt] = fl->brss;
00273 if (fl->num) {
00274 fte[Dpt] = fl->fte;
00275 fpth[Dpt] = fl->tt[fl->num];
00276 fpph[Dpt] = fl->pp[fl->num];
00277 } else {
00278 fte[Dpt] = 0.;
00279 failed++;
00280 }
00281
00282 }
00283 }
00284 free(fl); free(pt0);
00285
00286
00287 }
00288
00289
00290
00291
00292
00293
00294
00295 void fte_subearth(float *g, float *h, struct Gridsub *gridsub, float *fte_sub,
00296 float *fpph_sub, float *fpth_sub, float *brss_sub, int lmax,
00297 void (*integrator)(struct Point *, struct Point *, float *, float *, int, float, float,
00298 void (*)(float *, float *, struct Point *, float *, int, float)))
00299 {
00300 int np = gridsub->np, nt = gridsub->nt;
00301 int failed = 0;
00302
00303 struct Point *pt0 = (struct Point *)(malloc(sizeof(struct Point)));
00304 pt0->r = RS;
00305 struct FldLn *fl = (struct FldLn *)(malloc(sizeof(struct FldLn)));
00306 fl->num = 0;
00307 void (*Bfunc)(float *, float *, struct Point *, float *, int, float);
00308 Bfunc = pfss0;
00309
00310 for (int j = 0; j < 3; j++) {
00311 for (int i = 0; i < np; i++) {
00312 pt0->p = gridsub->ph[i];
00313 int Dpt = j * np + i;
00314 pt0->t = gridsub->th[Dpt];
00315
00316 trace(g, h, pt0, fl, lmax, integrator, Bfunc, RS, R0 + 0.001, 0.0);
00317
00318 brss_sub[Dpt] = fl->brss;
00319 if (fl->num) {
00320 fte_sub[Dpt] = fl->fte;
00321 fpth_sub[Dpt] = fl->tt[fl->num];
00322 fpph_sub[Dpt] = fl->pp[fl->num];
00323 } else {
00324 fte_sub[Dpt] = 0.;
00325 failed++;
00326 }
00327 }
00328 }
00329 free(fl); free(pt0);
00330
00331
00332 }
00333
00334
00335
00336
00337
00338
00339
00340 void propagate(double *t, float *v, float *b, double *t1AU, float *v1AU,
00341 float *b1AU, int np)
00342 {
00343 double step = PROPDIST * RSUN / NSTEP;
00344 double t1AU_old[np];
00345 float v1AU_old[np];
00346 float b1AU_old[np];
00347 for (int i = 0; i < np; i++) {
00348 t1AU[i] = t[i];
00349 v1AU[i] = v[i];
00350 b1AU[i] = b[i];
00351 }
00352 for (int rrr = 0; rrr < NSTEP; rrr++) {
00353 for (int i = 0; i < np; i++) {
00354 t1AU_old[i] = t1AU[i];
00355 v1AU_old[i] = v1AU[i];
00356 b1AU_old[i] = b1AU[i];
00357 }
00358
00359 for (int i = 0; i < np; i++) {
00360 t1AU[i] = t1AU_old[i] + (step / v1AU_old[i]);
00361 }
00362
00363 v1AU[0] = v1AU_old[0];
00364 b1AU[0] = b1AU_old[0];
00365 for (int qq = 1; qq < np; qq++) {
00366
00367 if (t1AU[qq] < t1AU[qq - 1]) {
00368 t1AU[qq] = t1AU[qq - 1] + 0.0001;
00369 double dtime = t1AU[qq] - t1AU_old[qq];
00370 float v0 = step / dtime;
00371 float v1 = v1AU_old[qq];
00372 v1AU[qq] = sqrt(2.0 / (1.0 / (v0 * v0) + 1.0 / (v1 * v1)));
00373 if (v1AU[qq] < v1AU[qq - 1])
00374 v1AU[qq] = v1AU[qq - 1];
00375 b1AU[qq] = (b1AU_old[qq] + b1AU_old[qq - 1]) / 2.0;
00376 } else {
00377 v1AU[qq] = v1AU_old[qq];
00378 b1AU[qq] = b1AU_old[qq];
00379 }
00380 }
00381 }
00382 float rar2 = (PROPDIST / RS) * (PROPDIST / RS);
00383 for (int i = 0; i < np; i++) b1AU[i] /= rar2;
00384 }
00385
00386
00387
00388 void interpol(float *y0, double *x0, int n0, float *y, double *x, int n)
00389 {
00390 if (n0 < 2) {
00391 printf("Not enough points in given series.\n");
00392 return;
00393 }
00394 int i1, i2;
00395 for (int i = 0; i < n; i++) {
00396 if (x[i] <= x0[0]) {
00397 i1 = 1; i2 = 0;
00398 } else if (x[i] > x0[n0 - 1]) {
00399 i1 = n0 - 2; i2 = n0 - 1;
00400 } else {
00401 for (int j = 0; j < n0 - 2; j++)
00402 if (x0[j] < x[i] && x0[j + 1] >= x[i]) {
00403 i1 = j; i2 = j + 1;
00404 break;
00405 }
00406 }
00407 y[i] = y0[i1] + (y0[i2] - y0[i1]) /
00408 (x0[i2] - x0[i1]) * (x[i] - x0[i1]);
00409 }
00410 }
00411
00412
00413
00414
00415
00416 void fixts1au(double *t, float *v, float *b, int n,
00417 double *date, float *speed, int *imf, int num)
00418 {
00419
00420 double t0[n];
00421 float v0[n], b0[n];
00422 for (int i = 0; i < n; i++) {
00423 t0[i] = 0; v0[i] = 0; b0[i] = 0;
00424 }
00425
00426 int n1 = 0, n2 = 1, n0 = 0;
00427 while (n1 < n) {
00428 while (n2 < n && fabs(t[n2] - t[n1]) < 90.)
00429 n2++;
00430 for (int i = n1; i < n2; i++) {
00431 t0[n0] += t[i]; v0[n0] += v[i]; b0[n0] += b[i];
00432 }
00433 int cc = n2 - n1;
00434 t0[n0] /= cc; v0[n0] /= cc; b0[n0] /= cc;
00435 n1 = n2; n2 = n1 + 1; n0++;
00436 }
00437
00438 float imf_t[num];
00439 interpol(v0, t0, n0, speed, date, num);
00440 interpol(b0, t0, n0, imf_t, date, num);
00441
00442 for (int i = 0; i < num; i++) {
00443 if (imf_t[i] > 0)
00444 imf[i] = -1;
00445 else if (imf_t[i] < 0)
00446 imf[i] = 1;
00447 else
00448 imf[i] = 0;
00449 }
00450 }
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461 void wsa_ss(float *g, float *h, int lmax, int rk4,
00462 struct Grid *grid, struct Gridsub *gridsub,
00463 float *v_g, float *brss_g, float *fte_g, float *fpph_g, float *fpth_g,
00464 float *v_s, float *brss_s, float *fte_s, float *fpph_s, float *fpth_s)
00465 {
00466 int np = grid->np, nt = grid->nt;
00467 int npnt = np * nt;
00468 int i, j, Dpt;
00469 int *fpmap;
00470 float ang;
00471
00472 void (*integrator)(struct Point *, struct Point *, float *, float *, int, float, float,
00473 void (*)(float *, float *, struct Point *, float *, int, float));
00474 if (rk4) integrator = RK4; else integrator = Euler;
00475
00476
00477 fte_global(g, h, grid, fte_g, fpph_g, fpth_g, brss_g, lmax, integrator);
00478
00479 fte_subearth(g, h, gridsub, fte_s, fpph_s, fpth_s, brss_s, lmax, integrator);
00480
00481 fpmap = (int *)(calloc(npnt, sizeof(int)));
00482 getbd(fpmap, fte_g, fpph_g, fpth_g, np, nt);
00483
00484 for (j = 0; j < nt; j++) {
00485 for (i = 0; i < np; i++) {
00486 Dpt = j * np + i;
00487 ang = ang_sep(fpph_g[Dpt], fpth_g[Dpt], fpmap, np, nt);
00488 v_g[Dpt] = wsa_f(fte_g[Dpt], ang);
00489 }
00490 }
00491
00492 for (j = 0; j < 3; j++) {
00493 for (i = 0; i < np; i++) {
00494 Dpt = j * np + i;
00495 ang = ang_sep(fpph_s[Dpt], fpth_s[Dpt], fpmap, np, nt);
00496 v_s[Dpt] = wsa_f(fte_s[Dpt], ang);
00497 }
00498 }
00499 free(fpmap);
00500 }
00501
00502
00503
00504
00505
00506
00507
00508
00509 void wsa_1AU(struct Gridsub *gridsub, float *v_s, float *brss_s, double *date,
00510 float *speed, int *imf, int tslen)
00511 {
00512 int i, j;
00513 int np = gridsub->np;
00514
00515
00516 double t0[np], t1au[np];
00517 float v0[np], b0[np], v1au[np], b1au[np];
00518 for (i = 0; i < np; i++)
00519 t0[i] = gridsub->time[np - 1 - i];
00520
00521
00522
00523
00524 float *speed_tmp;
00525 int *imf_tmp;
00526
00527 for (j = 0; j < 3; j++) {
00528 for (i = 0; i < np; i++) {
00529 v0[i] = v_s[j * np + np - 1 - i];
00530 b0[i] = brss_s[j * np + np - 1 - i];
00531 }
00532 propagate(t0, v0, b0, t1au, v1au, b1au, np);
00533 speed_tmp = speed + j * tslen; imf_tmp = imf + j * tslen;
00534 fixts1au(t1au, v1au, b1au, np, date, speed_tmp, imf_tmp, tslen);
00535 }
00536 }
00537