00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "cddefines.h"
00010 #include "physconst.h"
00011 #include "radius.h"
00012 #include "dense.h"
00013 #include "hyperfine.h"
00014 #include "magnetic.h"
00015 #include "hmi.h"
00016 #include "phycon.h"
00017 #include "geometry.h"
00018 #include "mean.h"
00019
00020
00021 void MeanInc(void)
00022 {
00023 long int ion,
00024 nelem;
00025 double dc,
00026 dcv;
00027
00028
00029
00030
00031 DEBUG_ENTRY( "MeanInc()" );
00032
00033 for( nelem=0; nelem < LIMELM; nelem++ )
00034 {
00035
00036 if( nelem==ipHYDROGEN )
00037 {
00038
00039
00040 ion = 2;
00041
00042
00043
00044
00045 dc = 2.*hmi.H2_total*radius.drad_x_fillfac;
00046 mean.xIonMeans[0][nelem][ion] += dc;
00047
00048 dc = 2.*hmi.H2_total*radius.drad_x_fillfac*dense.eden;
00049 mean.xIonEdenMeans[0][nelem][ion] += dc;
00050
00051 dcv = 2.*hmi.H2_total*radius.dVeff;
00052 mean.xIonMeans[1][nelem][ion] += dcv;
00053
00054 dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden;
00055 mean.xIonEdenMeans[1][nelem][ion] += dcv;
00056 }
00057 for( ion=0; ion < (nelem + 2); ion++ )
00058 {
00059
00060 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac;
00061 mean.xIonMeans[0][nelem][ion] += dc;
00062
00063
00064 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden;
00065 mean.xIonEdenMeans[0][nelem][ion] += dc;
00066
00067
00068
00069
00070 dcv = dense.xIonDense[nelem][ion]*radius.dVeff;
00071 mean.xIonMeans[1][nelem][ion] += dcv;
00072
00073
00074
00075
00076 dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden;
00077 mean.xIonEdenMeans[1][nelem][ion] += dcv;
00078 }
00079
00080
00081 dc = dense.gas_phase[nelem]*radius.drad_x_fillfac;
00082 mean.xIonMeansNorm[0][nelem] += dc;
00083 dc = dense.gas_phase[nelem]*radius.drad_x_fillfac*dense.eden;
00084 mean.xIonEdenMeansNorm[0][nelem] += dc;
00085
00086
00087
00088
00089
00090 dcv = dense.gas_phase[nelem]*radius.dVeff;
00091 mean.xIonMeansNorm[1][nelem] += dcv;
00092
00093 dcv = dense.gas_phase[nelem]*radius.dVeff*dense.eden;
00094 mean.xIonEdenMeansNorm[1][nelem] += dcv;
00095
00096
00097
00098
00099
00100
00101 if( nelem==ipHYDROGEN )
00102 {
00103 ion = 2;
00104
00105 dc = 2.*hmi.H2_total*radius.drad_x_fillfac;
00106 mean.TempMeans[0][nelem][ion] += dc*phycon.te;
00107 mean.TempMeansNorm[0][nelem][ion] += dc;
00108
00109
00110 dc = 2.*hmi.H2_total*radius.drad_x_fillfac*dense.eden;
00111 mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
00112 mean.TempEdenMeansNorm[0][nelem][ion] += dc;
00113
00114
00115
00116 dcv = 2.*hmi.H2_total*radius.dVeff;
00117 mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
00118 mean.TempMeansNorm[1][nelem][ion] += dcv;
00119
00120
00121 dcv = 2.*hmi.H2_total*radius.dVeff*dense.eden;
00122 mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te;
00123 mean.TempEdenMeansNorm[1][nelem][ion] += dcv;
00124 }
00125 for( ion=0; ion < (nelem + 2); ion++ )
00126 {
00127
00128 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac;
00129 mean.TempMeans[0][nelem][ion] += dc*phycon.te;
00130 mean.TempMeansNorm[0][nelem][ion] += dc;
00131
00132
00133 dc = dense.xIonDense[nelem][ion]*radius.drad_x_fillfac*dense.eden;
00134 mean.TempEdenMeans[0][nelem][ion] += dc*phycon.te;
00135 mean.TempEdenMeansNorm[0][nelem][ion] += dc;
00136
00137
00138
00139 dcv = dense.xIonDense[nelem][ion]*radius.dVeff;
00140 mean.TempMeans[1][nelem][ion] += dcv*phycon.te;
00141 mean.TempMeansNorm[1][nelem][ion] += dcv;
00142
00143
00144 dcv = dense.xIonDense[nelem][ion]*radius.dVeff*dense.eden;
00145 mean.TempEdenMeans[1][nelem][ion] += dcv*phycon.te;
00146 mean.TempEdenMeansNorm[1][nelem][ion] += dcv;
00147 }
00148 }
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 if( hyperfine.Tspin21cm > SMALLFLOAT )
00159 {
00160 dc = dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac/phycon.te;
00161 }
00162 else
00163 {
00164 dc = 0.;
00165 }
00166
00167 mean.B_HarMeanTempRadius[0] += sqrt(fabs(magnetic.pressure)*PI8) * dc;
00168
00169 mean.B_HarMeanTempRadius[1] += dc;
00170
00171
00172
00173 dc = dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac;
00174
00175
00176
00177
00178
00179 mean.HarMeanTempRadius[0] += dc;
00180 mean.HarMeanTempRadius[1] += dc/phycon.te;
00181
00182
00183 mean.HarMeanTempVolume[0] += dc*radius.r1r0sq;
00184 mean.HarMeanTempVolume[1] += dc/phycon.te*radius.r1r0sq;
00185
00186
00187 mean.H_21cm_spin_mean_radius[0] += dc;
00188 mean.H_21cm_spin_mean_radius[1] += dc / SDIV( hyperfine.Tspin21cm );
00189
00190
00191 dc = hmi.H2_total*radius.drad_x_fillfac;
00192 mean.H2MeanTempRadius[0] += dc*phycon.te;
00193 mean.H2MeanTempRadius[1] += dc;
00194
00195
00196 dcv = hmi.H2_total*radius.dVeff;
00197 mean.H2MeanTempVolume[0] += dc*phycon.te;
00198 mean.H2MeanTempVolume[1] += dc;
00199
00200
00201 dc = radius.drad_x_fillfac;
00202 mean.TempMeanRadius[0] += dc*phycon.te;
00203 mean.TempMeanRadius[1] += dc;
00204
00205 dcv = radius.dVeff;
00206 mean.TempMeanVolume[0] += dc*phycon.te;
00207 mean.TempMeanVolume[1] += dc;
00208
00209 DEBUG_EXIT( "MeanInc()" );
00210 return;
00211 }
00212
00213
00214 void MeanZero(void)
00215 {
00216 long int ion,
00217 j,
00218 nelem,
00219 limit;
00220 static bool lgFirst=true;
00221
00222 DEBUG_ENTRY( "MeanZero()" );
00223
00224
00225
00226
00227
00228
00229 if( lgFirst )
00230 {
00231 lgFirst = false;
00232
00233 mean.xIonMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
00234 mean.xIonEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
00235 mean.xIonMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) );
00236 mean.xIonEdenMeansNorm = (double **)MALLOC(sizeof(double *)*(unsigned)(2) );
00237 mean.TempMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
00238 mean.TempEdenMeans = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
00239 mean.TempMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
00240 mean.TempEdenMeansNorm = (double ***)MALLOC(sizeof(double **)*(unsigned)(2) );
00241 for( j=0; j<2; ++j )
00242 {
00243 mean.xIonMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
00244 mean.xIonEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
00245 mean.xIonMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) );
00246 mean.xIonEdenMeansNorm[j] = (double *)MALLOC(sizeof(double )*(unsigned)(LIMELM) );
00247 mean.TempMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
00248 mean.TempEdenMeans[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
00249 mean.TempMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
00250 mean.TempEdenMeansNorm[j] = (double **)MALLOC(sizeof(double *)*(unsigned)(LIMELM) );
00251 for( nelem=0; nelem<LIMELM; ++nelem )
00252 {
00253 limit = MAX2(3,nelem+2);
00254 mean.xIonMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
00255 mean.xIonEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
00256 mean.TempMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
00257 mean.TempEdenMeans[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
00258 mean.TempMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
00259 mean.TempEdenMeansNorm[j][nelem] = (double *)MALLOC(sizeof(double )*(unsigned)(limit) );
00260 }
00261 }
00262 }
00263
00264 for( j=0; j < 2; j++ )
00265 {
00266 for( nelem=0; nelem < LIMELM; nelem++ )
00267 {
00268 mean.xIonMeansNorm[j][nelem] = 0.;
00269 mean.xIonEdenMeansNorm[j][nelem] = 0.;
00270
00271 limit = MAX2(3,nelem+2);
00272
00273 for( ion=0; ion < limit; ion++ )
00274 {
00275 mean.TempMeansNorm[j][nelem][ion] = 0.;
00276 mean.TempEdenMeansNorm[j][nelem][ion] = 0.;
00277
00278 mean.xIonMeans[j][nelem][ion] = 0.;
00279 mean.xIonEdenMeans[j][nelem][ion] = 0.;
00280
00281 mean.TempMeans[j][nelem][ion] = 0.;
00282 mean.TempEdenMeans[j][nelem][ion] = 0.;
00283 }
00284 }
00285 }
00286
00287
00288 mean.HarMeanTempRadius[0] = 0.;
00289 mean.HarMeanTempRadius[1] = 0.;
00290
00291
00292 mean.HarMeanTempVolume[0] = 0.;
00293 mean.HarMeanTempVolume[1] = 0.;
00294
00295
00296 mean.H_21cm_spin_mean_radius[0] = 0.;
00297 mean.H_21cm_spin_mean_radius[1] = 0.;
00298
00299
00300 mean.H2MeanTempRadius[0] = 0.;
00301 mean.H2MeanTempRadius[1] = 0.;
00302
00303 mean.H2MeanTempVolume[0] = 0.;
00304 mean.H2MeanTempVolume[1] = 0.;
00305
00306 mean.TempMeanRadius[0] = 0.;
00307 mean.TempMeanRadius[1] = 0.;
00308 mean.TempMeanVolume[0] = 0.;
00309 mean.TempMeanVolume[1] = 0.;
00310
00311 mean.B_HarMeanTempRadius[0] = 0.;
00312 mean.B_HarMeanTempRadius[1] = 0.;
00313
00314 DEBUG_EXIT( "MeanZero()" );
00315 return;
00316 }
00317
00318
00319 void MeanIonRadius(
00320
00321 char chType,
00322
00323 long int nelem,
00324
00325 long int *n,
00326
00327
00328 float arlog[],
00329
00330 bool lgDensity )
00331 {
00332 long int ion,
00333 limit;
00334 double abund_radmean,
00335 normalize;
00336
00337 DEBUG_ENTRY( "MeanIonRadius()" );
00338
00339 ASSERT( chType=='i' || chType=='t' );
00340
00341
00342
00343
00344 limit = MAX2( 3, nelem+2 );
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354 if( !dense.lgElmtOn[nelem] )
00355 {
00356
00357
00358
00359 for( ion=0; ion < limit; ion++ )
00360 {
00361 arlog[ion] = -30.f;
00362 }
00363 *n = 0;
00364
00365 DEBUG_EXIT( "MeanIonRadius()" );
00366 return;
00367 }
00368
00369
00370
00371 *n = limit;
00372 while( *n > 0 && mean.xIonMeans[0][nelem][*n-1] <= 0. )
00373 {
00374 arlog[*n-1] = -30.f;
00375 --*n;
00376 }
00377
00378 if( chType=='i' && lgDensity)
00379 {
00380
00381 for( ion=0; ion < *n; ion++ )
00382 {
00383 if( mean.xIonEdenMeans[0][nelem][ion] > 0. )
00384 {
00385 abund_radmean = mean.xIonEdenMeans[0][nelem][ion];
00386 arlog[ion] = (float)log10(MAX2(1e-30,abund_radmean / mean.xIonEdenMeansNorm[0][nelem]));
00387 }
00388 else
00389 {
00390 arlog[ion] = -30.f;
00391 }
00392 }
00393 }
00394 else if( chType=='i' )
00395 {
00396
00397 for( ion=0; ion < *n; ion++ )
00398 {
00399 if( mean.xIonMeans[0][nelem][ion] > 0. )
00400 {
00401 abund_radmean = mean.xIonMeans[0][nelem][ion];
00402 arlog[ion] = (float)log10(MAX2(1e-30,abund_radmean/mean.xIonMeansNorm[0][nelem]));
00403 }
00404 else
00405 {
00406 arlog[ion] = -30.f;
00407 }
00408 }
00409 }
00410 else if( chType=='t' && lgDensity )
00411 {
00412
00413 for( ion=0; ion < *n; ion++ )
00414 {
00415 normalize = mean.TempEdenMeansNorm[0][nelem][ion];
00416 if( normalize > SMALLFLOAT )
00417 {
00418 abund_radmean = mean.TempEdenMeans[0][nelem][ion];
00419 arlog[ion] = (float)log10(MAX2(1e-30,abund_radmean/normalize));
00420 }
00421 else
00422 {
00423 arlog[ion] = -30.f;
00424 }
00425 }
00426 }
00427 else if( chType=='t' )
00428 {
00429
00430 for( ion=0; ion < *n; ion++ )
00431 {
00432 normalize = mean.TempMeansNorm[0][nelem][ion];
00433 if( normalize > SMALLFLOAT )
00434 {
00435 abund_radmean = mean.TempMeans[0][nelem][ion];
00436 arlog[ion] = (float)log10(MAX2(1e-30,abund_radmean/normalize));
00437 }
00438 else
00439 {
00440 arlog[ion] = -30.f;
00441 }
00442 }
00443 }
00444 else
00445 {
00446 fprintf(ioQQQ," MeanIonRadius called with insane job \n");
00447 }
00448
00449 DEBUG_EXIT( "MeanIonRadius()" );
00450 return;
00451 }
00452
00453
00454 void MeanIonVolume(
00455
00456 char chType,
00457
00458 long int nelem,
00459
00460 long int *n,
00461
00462 float arlog[],
00463
00464 bool lgDensity )
00465 {
00466 long int ion;
00467 double abund_volmean,
00468 normalize;
00469 long int limit;
00470
00471 DEBUG_ENTRY( "MeanIonVolume()" );
00472
00473 ASSERT( chType=='i' || chType=='t' );
00474
00475
00476
00477 limit = MAX2( 3, nelem+2 );
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488 if( !dense.lgElmtOn[nelem] )
00489 {
00490
00491 for( ion=0; ion <= limit; ion++ )
00492 {
00493 arlog[ion] = -30.f;
00494 }
00495 *n = 0;
00496
00497 DEBUG_EXIT( "MeanIonVolume()" );
00498 return;
00499 }
00500
00501
00502 *n = limit;
00503
00504
00505 while( *n > 0 && mean.xIonMeans[1][nelem][*n-1] <= 0. )
00506 {
00507 arlog[*n-1] = -30.f;
00508 --*n;
00509 }
00510
00511
00512 if( chType=='i' && lgDensity)
00513 {
00514
00515
00516
00517 for( ion=0; ion < *n; ion++ )
00518 {
00519 if( mean.xIonEdenMeans[1][nelem][ion] > 0. )
00520 {
00521 abund_volmean = mean.xIonEdenMeans[1][nelem][ion];
00522 arlog[ion] = (float)log10(MAX2(1e-30,abund_volmean) / mean.xIonEdenMeansNorm[1][nelem]);
00523 }
00524 else
00525 {
00526 arlog[ion] = -30.f;
00527 }
00528 }
00529 }
00530 else if( chType=='i' )
00531 {
00532
00533
00534
00535 for( ion=0; ion < *n; ion++ )
00536 {
00537 if( mean.xIonMeans[1][nelem][ion] > 0. )
00538 {
00539 abund_volmean = mean.xIonMeans[1][nelem][ion];
00540 arlog[ion] = (float)log10(MAX2(1e-30,abund_volmean) / mean.xIonMeansNorm[1][nelem]);
00541 }
00542 else
00543 {
00544 arlog[ion] = -30.f;
00545 }
00546 }
00547 }
00548 else if( chType=='t' && lgDensity )
00549 {
00550
00551
00552
00553 for( ion=0; ion < *n; ion++ )
00554 {
00555 normalize = mean.TempEdenMeansNorm[1][nelem][ion];
00556 if( normalize > SMALLFLOAT )
00557 {
00558 abund_volmean = mean.TempEdenMeans[1][nelem][ion];
00559 arlog[ion] = (float)log10(MAX2(1e-30,abund_volmean)/normalize);
00560 }
00561 else
00562 {
00563 arlog[ion] = -30.f;
00564 }
00565 }
00566 }
00567 else if( chType=='t' )
00568 {
00569
00570
00571
00572 for( ion=0; ion < *n; ion++ )
00573 {
00574 normalize = mean.TempMeansNorm[1][nelem][ion];
00575 if( normalize > SMALLFLOAT )
00576 {
00577 abund_volmean = mean.TempMeans[1][nelem][ion];
00578 arlog[ion] = (float)log10(MAX2(1e-30,abund_volmean)/normalize);
00579 }
00580 else
00581 {
00582 arlog[ion] = -30.f;
00583 }
00584 }
00585 }
00586 else
00587 {
00588 fprintf(ioQQQ," MeanIonVolume called with insane job\n");
00589 }
00590
00591 DEBUG_EXIT( "MeanIonVolume()" );
00592 return;
00593 }
00594
00595
00596
00597 void aver(
00598
00599 const char *chWhat,
00600
00601 double quan,
00602
00603 double weight,
00604
00605 const char *chLabl)
00606 {
00607
00608 # define NAVER 20
00609 static char chLabavr[NAVER][11];
00610 long int i,
00611 ioff,
00612 j;
00613 static long int n;
00614 double raver[NAVER]={0.},
00615 vaver[NAVER]={0.};
00616 static double aversv[4][NAVER];
00617
00618 DEBUG_ENTRY( "aver()" );
00619
00620 if( strcmp(chWhat,"zero") == 0 )
00621 {
00622
00623 for( i=0; i < NAVER; i++ )
00624 {
00625 for( j=0; j < 4; j++ )
00626 {
00627 aversv[j][i] = 0.;
00628 }
00629 }
00630 n = 0;
00631 }
00632 else if( strcmp(chWhat,"zone") == 0 )
00633 {
00634 n = 0;
00635 }
00636 else if( strcmp(chWhat,"doit") == 0 )
00637 {
00638
00639 if( n >= NAVER )
00640 {
00641 fprintf( ioQQQ, " Too many values entered into AVER, increase NAVER\n" );
00642 puts( "[Stop in aver]" );
00643 cdEXIT(EXIT_FAILURE);
00644 }
00645 aversv[0][n] += weight*quan*radius.drad_x_fillfac;
00646 aversv[1][n] += weight*radius.drad_x_fillfac;
00647 aversv[2][n] += weight*quan*radius.dVeff;
00648 aversv[3][n] += weight*radius.dVeff;
00649
00650
00651 strcpy( chLabavr[n], chLabl );
00652 n += 1;
00653 }
00654 else if( strcmp(chWhat,"prin") == 0 )
00655 {
00656
00657 ioff = n*10/2 + 1;
00658
00659
00660 fprintf( ioQQQ,"\n");
00661 for( i=0; i < ioff; i++ )
00662 {
00663
00664 fprintf( ioQQQ, " " );
00665 }
00666
00667
00668 fprintf( ioQQQ, "Averaged Quantities\n " );
00669
00670 fprintf( ioQQQ, " " );
00671 for( i=0; i < n; i++ )
00672 {
00673 fprintf( ioQQQ, "%10.10s", chLabavr[i] );
00674 }
00675 fprintf( ioQQQ, " \n" );
00676
00677 for( i=0; i < n; i++ )
00678 {
00679 if( aversv[1][i] > 0. )
00680 {
00681 raver[i] = aversv[0][i]/aversv[1][i];
00682 }
00683 else
00684 {
00685 raver[i] = 0.;
00686 }
00687 if( aversv[3][i] > 0. )
00688 {
00689 vaver[i] = aversv[2][i]/aversv[3][i];
00690 }
00691 else
00692 {
00693 vaver[i] = 0.;
00694 }
00695 }
00696
00697 fprintf( ioQQQ, " Radius:" );
00698 for( i=0; i < n; i++ )
00699 {
00700
00701 fprintf( ioQQQ, " ");
00702 fprintf(ioQQQ,PrintEfmt("%9.2e", raver[i] ) );
00703 }
00704 fprintf( ioQQQ, "\n" );
00705
00706
00707 if( geometry.lgSphere )
00708 {
00709 fprintf( ioQQQ, " Volume:" );
00710 for( i=0; i < n; i++ )
00711 {
00712
00713 fprintf( ioQQQ, " ");
00714 fprintf(ioQQQ,PrintEfmt("%9.2e", vaver[i] ) );
00715 }
00716 fprintf( ioQQQ, "\n" );
00717 }
00718 }
00719 else
00720 {
00721 fprintf( ioQQQ, " AVER does not understand chWhat=%4.4s\n",
00722 chWhat );
00723 ShowMe();
00724 puts( "[Stop in aver]" );
00725 cdEXIT(EXIT_FAILURE);
00726 }
00727
00728 DEBUG_EXIT( "aver()" );
00729 return;
00730 }
00731