00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "cddefines.h"
00012 #include "taulines.h"
00013 #include "opacity.h"
00014 #include "dense.h"
00015 #include "iso.h"
00016 #include "hmi.h"
00017 #include "h2.h"
00018 #include "rfield.h"
00019 #include "conv.h"
00020 #include "rt.h"
00021 #include "mole_co_atom.h"
00022 #include "atomfeii.h"
00023 #include "heavy.h"
00024 #include "he.h"
00025 #include "trace.h"
00026
00027
00028 static int nOTS_Line_type = -1;
00029 static int nOTS1=-1 , nOTS2=-1;
00030
00031 static void RT_OTS_AddCont(
00032
00033 float ots,
00034
00035 long int ip);
00036
00037
00038
00039 void RT_OTS(void)
00040 {
00041 long int
00042 ipla,
00043 ipISO ,
00044 nelem,
00045 n;
00046 float
00047 difflya,
00048 esc,
00049 ots;
00050
00051
00052
00053 # define BOWEN 0.5f
00054 long int ipHi,
00055 ipLo;
00056
00057 double bwnfac;
00058 double ots660;
00059 float cont_phot_destroyed;
00060 double save_lya_dest,
00061 save_he2lya_dest;
00062
00063 double save_he2rec_dest;
00064
00065
00066
00067
00068
00069 DEBUG_ENTRY( "RT_OTS()" );
00070
00071
00072
00073
00074
00075
00076 nOTS_Line_type = 0;
00077 nelem = ipHELIUM;
00078 if( dense.lgElmtOn[nelem] )
00079 {
00080
00081
00082 bwnfac = BOWEN * MAX2(0.f,1.f- EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pesc -
00083 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pelec_esc );
00084
00085
00086
00087 ots660 = EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Aul*
00088 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].PopHi*dense.xIonDense[nelem][nelem+1]*
00089
00090 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest *BOWEN*2.0;
00091
00092
00093 if( ots660 > SMALLFLOAT )
00094 RT_OTS_AddLine(ots660 , he.ip660 );
00095
00096
00097
00098 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest *= (float)bwnfac;
00099 ASSERT( EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].Pdest >= 0. );
00100 {
00101
00102
00103
00104 enum {DEBUG_LOC=false};
00105
00106 if( DEBUG_LOC )
00107 {
00108 fprintf(ioQQQ,"DEBUG HeII Bowen nzone %li bwnfac:%.2e bwnfac esc:%.2e ots660 %.2e\n",
00109 nzone,
00110 bwnfac ,
00111 bwnfac/BOWEN ,
00112 ots660 );
00113 }
00114 }
00115 }
00116
00117 else
00118 {
00119 bwnfac = 1.;
00120 }
00121
00122
00123 save_lya_dest = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest;
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest *= rfield.lgLyaOTS;
00137
00138
00139 save_he2lya_dest = EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest;
00140 EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest *= rfield.lgHeIIOTS;
00141
00142
00143 save_he2rec_dest = iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad];
00144 iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] *= rfield.lgHeIIOTS;
00145
00146 nOTS_Line_type = 1;
00147
00148
00149
00150
00151 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00152 {
00153 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00154 {
00155 nOTS1 = ipISO;
00156 nOTS2 = nelem;
00157
00161 if( (dense.IonHigh[nelem] >= nelem+1-ipISO ) )
00162 {
00163
00164
00165
00166
00167 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
00168 {
00169 for( ipLo=0; ipLo < ipHi; ipLo++ )
00170 {
00171
00172
00173
00174 if( EmisLines[ipISO][nelem][ipHi][ipLo].ipCont< 1 ||
00175 EmisLines[ipISO][nelem][ipHi][ipLo].Pdest<= DEST0 )
00176 continue;
00177
00178
00179 EmisLines[ipISO][nelem][ipHi][ipLo].ots =
00180 EmisLines[ipISO][nelem][ipHi][ipLo].Aul*
00181 EmisLines[ipISO][nelem][ipHi][ipLo].PopHi*dense.xIonDense[nelem][nelem+1-ipISO]*
00182 EmisLines[ipISO][nelem][ipHi][ipLo].Pdest;
00183
00184 ASSERT( EmisLines[ipISO][nelem][ipHi][ipLo].ots >= 0. );
00185
00186
00187
00188
00189
00190
00191
00192 if( EmisLines[ipISO][nelem][ipHi][ipLo].ots > SMALLFLOAT )
00193 RT_OTS_AddLine(EmisLines[ipISO][nelem][ipHi][ipLo].ots,
00194 EmisLines[ipISO][nelem][ipHi][ipLo].ipCont );
00195 }
00196 }
00197 {
00198
00199
00200
00201 enum {DEBUG_LOC=false};
00202
00203 if( DEBUG_LOC )
00204 {
00205 long ip;
00206 if( ipISO==0 && nelem==0 && nzone>500 )
00207 {
00208 ipHi = 2;
00209 ipLo = 0;
00210 ip = EmisLines[ipISO][nelem][ipHi][ipLo].ipCont;
00211 fprintf(ioQQQ,"DEBUG hlyaots\t%.2f\tEdenTrue\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00212 fnzone,
00213 dense.EdenTrue ,
00214 EmisLines[ipISO][nelem][ipHi][ipLo].ots,
00215 opac.opacity_abs[ip-1],
00216 EmisLines[ipISO][nelem][ipHi][ipLo].PopHi*dense.xIonDense[nelem][nelem+1-ipISO],
00217 EmisLines[ipISO][nelem][ipHi][ipLo].PopHi,
00218 dense.xIonDense[nelem][nelem+1-ipISO],
00219 EmisLines[ipISO][nelem][ipHi][ipLo].Pdest,
00220 rfield.otslinNew[ip-1],
00221 rfield.otslin[ip-1]);
00222 }
00223 }
00224 }
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
00235 {
00236 cont_phot_destroyed = (float)(iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
00237 (1. - iso.RadRecomb[ipISO][nelem][n][ipRecEsc])*dense.eden*
00238 dense.xIonDense[nelem][nelem+1]);
00239 ASSERT( cont_phot_destroyed >= 0. );
00240
00241
00242 RT_OTS_AddCont(cont_phot_destroyed,iso.ipIsoLevNIonCon[ipISO][nelem][n]);
00243
00244 {
00245
00246 enum {DEBUG_LOC=false};
00247
00248 if( DEBUG_LOC && nzone > 400 && nelem==0 && n==2 )
00249 {
00250 long ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1;
00251 fprintf(ioQQQ,"hotsdebugg\t%.3f\t%li\th con ots\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00252 fnzone,
00253 n ,
00254 EmisLines[ipISO][ipHYDROGEN][2][0].PopHi*dense.xIonDense[ipHYDROGEN][1],
00255 cont_phot_destroyed,
00256 cont_phot_destroyed/opac.opacity_abs[ip],
00257 rfield.otsconNew[ip] ,
00258 rfield.otscon[ip] ,
00259 opac.opacity_abs[ip] ,
00260 hmi.Hmolec[ipMHm] ,
00261 hmi.HMinus_photo_rate);
00262 }
00263 }
00264 }
00265 }
00266 }
00267 }
00268
00269 {
00270
00271 enum {DEBUG_LOC=false};
00272
00273 if( DEBUG_LOC )
00274 {
00275 nelem = 0;
00276 fprintf(ioQQQ,"hotsdebugg %li \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00277 nzone,
00278 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][0]-1],
00279 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][1]-1],
00280 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][3]-1],
00281 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][4]-1],
00282 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][5]-1],
00283 rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1],
00284 opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1]);
00285 }
00286 }
00287
00288
00289 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest = (float)save_lya_dest;
00290 EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest = (float)save_he2lya_dest;
00291 iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] = save_he2rec_dest;
00292
00293 nelem = ipHELIUM;
00294 if( dense.lgElmtOn[nelem] && bwnfac > 0. )
00295 {
00296
00297 EmisLines[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Pdest /= (float)bwnfac;
00298 }
00299
00300 if( trace.lgTrace )
00301 {
00302 fprintf(ioQQQ," RT_OTS Pdest %.2e ots rate %.2e in otslinNew %.2e con opac %.2e\n",
00303 EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest , EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ots ,
00304 rfield.otslinNew[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ,
00305 opac.opacity_abs[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1]
00306 );
00307 }
00308
00309 nOTS_Line_type = 2;
00310
00311 for( nelem=NISO; nelem < LIMELM; nelem++ )
00312 {
00313 long int ion;
00314
00315
00316 for( ion=0; ion < nelem+1-NISO; ion++ )
00317 {
00318 if( dense.xIonDense[nelem][ion+1] > 0. )
00319 {
00320 nOTS1 = nelem;
00321 nOTS2 = ion;
00322
00323 ipla = Heavy.ipLyHeavy[nelem][ion];
00324 esc = opac.ExpmTau[ipla-1];
00325
00326 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
00327
00328
00329 ots = difflya*MAX2(0.f,1.f-esc);
00330
00331
00332 ASSERT( ots >= 0.);
00333
00334
00335
00336
00337 if( ots > SMALLFLOAT )
00338 RT_OTS_AddLine(ots,ipla);
00339
00340
00341 ipla = Heavy.ipBalHeavy[nelem][ion];
00342 esc = opac.ExpmTau[ipla-1];
00343
00344 difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
00345
00346
00347 ots = difflya*MAX2(0.f,1.f-esc);
00348 ASSERT( ots >= 0.);
00349
00350
00351 if( ots > SMALLFLOAT )
00352 RT_OTS_AddLine(ots,ipla);
00353 }
00354 }
00355 }
00356
00357 nOTS_Line_type = 3;
00358
00359 FeII_OTS();
00360
00361 nOTS_Line_type = 4;
00362
00363 # if 0
00364 {
00365 # include "lines_service.h"
00366 if( nzone>290 ) DumpLine( &TauLines[74] );
00367 }
00368 # endif
00369
00370
00371 for( nOTS1=1; nOTS1 < nLevel1; nOTS1++ )
00372 {
00373
00374 TauLines[nOTS1].ots = TauLines[nOTS1].PopHi * TauLines[nOTS1].Aul * TauLines[nOTS1].Pdest;
00375
00376
00377
00378 if( TauLines[nOTS1].ots > SMALLFLOAT )
00379 RT_OTS_AddLine( TauLines[nOTS1].ots , TauLines[nOTS1].ipCont);
00380 }
00381
00382 nOTS_Line_type = 5;
00383
00384 for( nOTS1=0; nOTS1 < nWindLine; nOTS1++ )
00385 {
00386 if( TauLine2[nOTS1].IonStg < TauLine2[nOTS1].nelem+1-NISO )
00387 {
00388 TauLine2[nOTS1].ots = TauLine2[nOTS1].PopHi * TauLine2[nOTS1].Aul * TauLine2[nOTS1].Pdest;
00389 if( TauLine2[nOTS1].ots > SMALLFLOAT )
00390 RT_OTS_AddLine( TauLine2[nOTS1].ots , TauLine2[nOTS1].ipCont);
00391 }
00392 }
00393
00394
00395
00396
00397
00398
00399 nOTS_Line_type = 6;
00400 CO_OTS();
00401
00402
00403
00404 nOTS_Line_type = 7;
00405 H2_RT_OTS();
00406
00407
00408 DEBUG_EXIT( "RT_OTS()" );
00409
00410 return;
00411 }
00412
00413
00414
00415 void RT_OTS_AddLine(double ots,
00416
00417 long int ip )
00418 {
00419
00420 DEBUG_ENTRY( "RT_OTS_AddLine()" );
00421
00422
00423
00424
00425
00426 if( ip==0 || ip > rfield.nflux )
00427 {
00428 DEBUG_EXIT( "RT_OTS_AddLine()" );
00429 return;
00430 }
00431
00432
00433 ASSERT( ots >= 0. );
00434
00435 ASSERT( ip > 0 );
00436
00437
00438
00439
00440 if( opac.opacity_abs[ip-1] > 0. )
00441 {
00442 rfield.otslinNew[ip-1] += (float)(ots/opac.opacity_abs[ip-1]);
00443 }
00444
00445 {
00446
00447 enum {DEBUG_LOC=false};
00448
00449 if( DEBUG_LOC && ip== 2363 )
00450 {
00451 fprintf(ioQQQ,"DEBUG ots, opc, otsr %.3e\t%.3e\t%.3e\t",
00452 ots ,
00453 opac.opacity_abs[ip-1],
00454 ots/opac.opacity_abs[ip-1] );
00455 fprintf(ioQQQ,"iteration %li type %i %i %i \n",
00456 iteration,
00457 nOTS_Line_type,
00458 nOTS1,nOTS2 );
00459 }
00460 }
00461
00462 DEBUG_EXIT( "RT_OTS_AddLine()" );
00463 return;
00464 }
00465
00466
00467
00468
00469 static void RT_OTS_AddCont(
00470
00471 float ots,
00472
00473 long int ip)
00474 {
00475
00476 DEBUG_ENTRY( "RT_OTS_AddCont()" );
00477
00478
00479
00480
00481
00482
00483
00484 if( ip > rfield.nflux )
00485 {
00486 DEBUG_EXIT( "RT_OTS_AddCont()" );
00487 return;
00488 }
00489
00490 ASSERT( ip> 0 );
00491
00492 ASSERT( ots >= 0. );
00493
00494 ASSERT( ip <= rfield.nupper );
00495
00496
00497
00498
00499 if( opac.opacity_abs[ip-1] > 0. )
00500 {
00501 rfield.otsconNew[ip-1] += (float)(ots/opac.opacity_abs[ip-1]);
00502
00503
00504
00505 }
00506
00507 DEBUG_EXIT( "RT_OTS_AddCont()" );
00508 return;
00509 }
00510
00511 #if 0
00512
00513 static bool lgLyaOTS_oscil( float ots_new )
00514 {
00515 static float old , older;
00516 bool lgReturn = false;
00517
00518 if( !conv.nTotalIoniz || conv.lgSearch )
00519 {
00520
00521 old = 0.;
00522 older = 0.;
00523 }
00524
00525
00526 if( conv.nPres2Ioniz>1 && (ots_new-old) * (old-older) < 0.)
00527 lgReturn = true;
00528
00529 older = old;
00530 old = ots_new;
00531 return( lgReturn );
00532 }
00533
00534
00535 if( lgLyaOTS_oscil( rfield.otslinNew[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ) )
00536 {
00537 fprintf(ioQQQ,"DEBUG lya ots %.2f %.3e\n",
00538 fnzone,
00539 rfield.otslinNew[EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
00540 }
00541 #endif
00542
00543
00544
00545 void RT_OTS_Update(
00546
00547 double *SumOTS ,
00548
00549
00550
00551 double BigFrac)
00552 {
00553 long int i;
00554
00555 static bool lgNeedSpace=true;
00556 static double *old_OTS_line_x_opac , *old_OTS_cont_x_opac;
00557 float FacBig , FacSml;
00558
00559 DEBUG_ENTRY( "RT_OTS_Update()" );
00560
00561 if( BigFrac <= 0. )
00562 BigFrac = 10.;
00563 FacBig = (float)(1. + BigFrac);
00564 FacSml = 1.f/FacBig;
00565
00566
00567 if( lgNeedSpace )
00568 {
00569 old_OTS_line_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) );
00570 old_OTS_cont_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) );
00571 }
00572 lgNeedSpace = false;
00573
00574 *SumOTS = 0.;
00575
00576
00577 if( rfield.lgKillOTSLine )
00578 {
00579 for( i=0; i < rfield.nflux; i++ )
00580 {
00581 rfield.otslinNew[i] = 0.;
00582 }
00583 }
00584
00585
00586 if( nzone==0 )
00587 {
00588 for( i=0; i < rfield.nflux; i++ )
00589 {
00590 old_OTS_line_x_opac[i] = rfield.otslinNew[i] * opac.opacity_abs[i];
00591 old_OTS_cont_x_opac[i] = rfield.otsconNew[i] * opac.opacity_abs[i];
00592 }
00593 }
00594
00595
00596
00597 *SumOTS = 0.;
00598
00599 for( i=0; i < rfield.nflux; ++i )
00600 {
00601 double OTS_line_x_opac = rfield.otslinNew[i] * opac.opacity_abs[i];
00602 double OTS_cont_x_opac = rfield.otsconNew[i] * opac.opacity_abs[i];
00603 double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );
00604
00605
00606 if( OTS_line_x_opac > old_OTS_line_x_opac[i]*FacBig )
00607 {
00608 rfield.otslin[i] = rfield.otslin[i]*FacBig;
00609 }
00610 else if( OTS_line_x_opac < old_OTS_line_x_opac[i]*FacSml )
00611 {
00612 rfield.otslin[i] = rfield.otslin[i]*FacSml;
00613 }
00614 else
00615 {
00616 rfield.otslin[i] = rfield.otslinNew[i];
00617 }
00618
00619
00620 if( OTS_cont_x_opac > old_OTS_cont_x_opac[i]*FacBig )
00621 {
00622 rfield.otscon[i] = rfield.otscon[i]*FacBig;
00623 }
00624 else if( OTS_cont_x_opac < old_OTS_cont_x_opac[i]*FacSml )
00625 {
00626 rfield.otscon[i] = rfield.otscon[i]*FacSml;
00627 }
00628 else
00629 {
00630 rfield.otscon[i] = rfield.otsconNew[i];
00631 }
00632
00633
00634 rfield.otsconNew[i] = 0.;
00635
00636
00637 rfield.otslinNew[i] = 0.;
00638
00639
00640 old_OTS_line_x_opac[i] = rfield.otslin[i] * opac.opacity_abs[i];
00641 old_OTS_cont_x_opac[i] = rfield.otscon[i] * opac.opacity_abs[i];
00642
00643
00644
00645 rfield.ConOTS_local_OTS_rate[i] = (float)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655 *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
00656
00657 rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+
00658 rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] +
00659 rfield.ConOTS_local_OTS_rate[i];
00660
00661 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
00662 rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
00663
00664 }
00665
00666
00667
00668
00669 rfield.flux_accum[rfield.nflux-1] = 0;
00670 for( i=1; i < rfield.nflux; i++ )
00671 {
00672 rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] +
00673 rfield.SummedCon[rfield.nflux-i-1];
00674 }
00675
00676 ASSERT( rfield.flux_accum[0]> 0. );
00677
00678
00679
00680 ASSERT(rfield.ipPlasma>0 );
00681
00682
00683 for( i=0; i < rfield.ipPlasma-1; i++ )
00684 {
00685 rfield.otscon[i] = 0.;
00686 rfield.ConOTS_local_OTS_rate[i] = 0.;
00687 rfield.ConOTS_local_photons[i] = 0.;
00688 rfield.otslin[i] = 0.;
00689 rfield.SummedDif[i] = 0.;
00690 rfield.OccNumbBremsCont[i] = 0.;
00691 rfield.SummedCon[i] = 0.;
00692 rfield.SummedOcc[i] = 0.;
00693 rfield.ConInterOut[i] = 0.;
00694 }
00695
00696
00697 if( rfield.ipEnergyBremsThin > 0 )
00698 {
00699 for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
00700 {
00701
00702
00703
00704
00705 float factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);
00706
00707
00708 rfield.OccNumbBremsCont[i] = (float)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
00709 }
00710 }
00711
00712 DEBUG_EXIT( "RT_OTS_Update()" );
00713 return;
00714
00715 # if 0
00716 OTSOld = 0.;
00717 OTSNew = 0.;
00718 BigOTSNew = 0.;
00719
00720 for( i=0; i < rfield.nflux; i++ )
00721 {
00722 OTSOld += rfield.otscon[i] + rfield.otslin[i];
00723 OTSNew += rfield.otsconNew[i] + rfield.otslinNew[i];
00724 if( BigOTSNew < rfield.otsconNew[i] + rfield.otslinNew[i] )
00725 {
00726 BigOTSNew = rfield.otsconNew[i] + rfield.otslinNew[i];
00727 }
00728 }
00729
00730
00731
00732 if( BigFrac == 0. || conv.lgSearch || OTSOld < SMALLFLOAT )
00733 {
00734
00735
00736 FracOld = 0.;
00737 }
00738 else
00739 {
00740 if( OTSNew > OTSOld )
00741 {
00742 chng = fabs(1. - OTSOld/SDIV(OTSNew) );
00743 }
00744 else
00745 {
00746 chng = fabs(1. - OTSNew/SDIV(OTSOld) );
00747 }
00748
00749 if( chng < BigFrac )
00750 {
00751
00752 FracOld = 0.;
00753 }
00754 else
00755 {
00756
00757 FracOld = (1. - BigFrac / chng);
00758 ASSERT( FracOld >= 0. );
00759 FracOld = MIN2( 0.25 , FracOld );
00760 }
00761 }
00762
00763
00764 FracNew = 1. - FracOld;
00765 fprintf(ioQQQ," DEBUG FracNew\t%.2e\n", FracNew);
00766
00767
00768 changeOTS = 0.;
00769
00770
00771 for( i=0; i < rfield.nflux; i++ )
00772 {
00773
00774 double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );
00775
00776
00777 if( fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i])> changeOTS)
00778 {
00779 changeOTS =
00780 fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i]);
00781 }
00782
00783
00784
00785
00786 if( BigFrac > 0. && conv.nTotalIoniz > 0 )
00787 {
00788
00789
00790
00791
00792
00793 rfield.otscon[i] = (float)((rfield.otscon[i]*FracOld +rfield.otsconNew[i]*FracNew));
00794
00795
00796
00797
00798 rfield.ConOTS_local_OTS_rate[i] = (float)(rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);
00799
00800
00801 rfield.otslin[i] = (float)((rfield.otslin[i]*FracOld + rfield.otslinNew[i]*FracNew));
00802 }
00803
00804
00805 rfield.otsconNew[i] = 0.;
00806
00807
00808 rfield.otslinNew[i] = 0.;
00809
00810
00811 *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
00812 {
00813
00814
00815 enum {DEBUG_LOC=false};
00816
00817 if( DEBUG_LOC && (nzone==378) )
00818 {
00819 if( conv.nPres2Ioniz > 3 )
00820 exit(1);
00821 fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
00822 conv.nPres2Ioniz,
00823 rfield.anu[i],
00824 opac.opacity_abs[i],
00825 rfield.otscon[i],
00826 rfield.otslin[i]);
00827 }
00828 }
00829
00830
00831
00832 rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+
00833 rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] +
00834 rfield.ConOTS_local_OTS_rate[i];
00835
00836 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
00837 rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
00838 }
00839
00840 ASSERT(*SumOTS >= 0. );
00841
00842
00843
00844 rfield.flux_accum[rfield.nflux-1] = 0;
00845 for( i=1; i < rfield.nflux; i++ )
00846 {
00847 rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] +
00848 rfield.SummedCon[rfield.nflux-i-1];
00849 }
00850
00851 ASSERT( rfield.flux_accum[0]> 0. );
00852
00853
00854
00855 ASSERT(rfield.ipPlasma>0 );
00856
00857
00858 for( i=0; i < rfield.ipPlasma-1; i++ )
00859 {
00860 rfield.otscon[i] = 0.;
00861 rfield.ConOTS_local_OTS_rate[i] = 0.;
00862 rfield.ConOTS_local_photons[i] = 0.;
00863 rfield.otslin[i] = 0.;
00864 rfield.SummedDif[i] = 0.;
00865 rfield.OccNumbBremsCont[i] = 0.;
00866 rfield.SummedCon[i] = 0.;
00867 rfield.SummedOcc[i] = 0.;
00868 rfield.ConInterOut[i] = 0.;
00869 }
00870
00871
00872 if( rfield.ipEnergyBremsThin > 0 )
00873 {
00874 for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
00875 {
00876
00877
00878
00879
00880 float factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);
00881
00882
00883 rfield.OccNumbBremsCont[i] = (float)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
00884 }
00885 }
00886
00887 {
00888
00889
00890 enum {DEBUG_LOC=false};
00891
00892 if( DEBUG_LOC && (nzone>0) )
00893 {
00894 double BigOTS;
00895 int ipOTS=0;
00896 BigOTS = -1.;
00897
00898
00899 for( i=0; i < rfield.nflux; i++ )
00900 {
00901 if( (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i] > BigOTS )
00902 {
00903 BigOTS = (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
00904 ipOTS = (int)i;
00905 }
00906 }
00907 fprintf(ioQQQ,
00908 " sumcontinuunots zone %li SumOTS %.2e %.2eRyd opc:%.2e OTS%.2e lin%.2e con:%.2e %s %s cell%i FracNew %.2f \n",
00909 nzone,
00910 *SumOTS,
00911 rfield.anu[ipOTS] ,
00912 opac.opacity_abs[ipOTS],
00913 BigOTS ,
00914 rfield.otslin[ipOTS],
00915 rfield.otscon[ipOTS] ,
00916
00917 rfield.chLineLabel[ipOTS] ,
00918
00919 rfield.chContLabel[ipOTS] ,
00920 ipOTS,
00921 FracNew);
00922 }
00923 }
00924
00925 DEBUG_EXIT( "RT_OTS_Update()" );
00926 return;
00927 # endif
00928 }
00929
00930
00931
00932
00933 void RT_OTS_Zero( void )
00934 {
00935 long int i;
00936
00937 DEBUG_ENTRY( "RT_OTS_Zero()" );
00938
00939
00940
00941 for( i=0; i <= rfield.nflux; i++ )
00942 {
00943 rfield.SummedDif[i] = 0.;
00944
00945 rfield.otscon[i] = 0.;
00946 rfield.otslin[i] = 0.;
00947 rfield.otslinNew[i] = 0.;
00948 rfield.otsconNew[i] = 0.;
00949
00950 rfield.trans_coef_zone[i] = 1.f;
00951 rfield.trans_coef_total[i] = 1.f;
00952
00953 rfield.ConInterOut[i] = 0.;
00954 rfield.outlin[i] = 0.;
00955 rfield.outlin_noplot[i] = 0.;
00956 rfield.SummedDif[i] = 0.;
00957
00958 rfield.SummedCon[i] = rfield.flux[i];
00959 rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
00960 rfield.ConOTS_local_photons[i] = 0.;
00961 rfield.ConOTS_local_OTS_rate[i] = 0.;
00962 }
00963
00964 DEBUG_EXIT( "RT_OTS_Zero()" );
00965 return;
00966 }
00967
00968
00969
00970
00971 void RT_OTS_ChkSum(long int ipPnt)
00972 {
00973 static long int nInsane=0;
00974 long int i;
00975 double phisig;
00976 # define LIM_INSAME_PRT 30
00977
00978 DEBUG_ENTRY( "RT_OTS_ChkSum()" );
00979
00980
00981
00982
00983 for( i=rfield.ipEnergyBremsThin; i < rfield.nflux; i++ )
00984 {
00985 phisig = rfield.otscon[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly +
00986 rfield.outlin[i]+
00987 rfield.outlin_noplot[i]+
00988 rfield.ConOTS_local_OTS_rate[i];
00989
00990
00991 if( phisig > 0. && rfield.SummedDif[i]> 0.)
00992 {
00993 if( fabs(rfield.SummedDif[i]/phisig-1.) > 1e-3 )
00994 {
00995 ++nInsane;
00996
00997
00998 if( nInsane < LIM_INSAME_PRT )
00999 {
01000 fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum insane SummedDif at energy %.5e error= %.2e i=%4ld\n",
01001 rfield.anu[i], rfield.SummedDif[i]/phisig - 1., i );
01002 fprintf( ioQQQ, " SummedDif, sum are%11.4e%11.4e\n",
01003 rfield.SummedDif[i], phisig );
01004 fprintf( ioQQQ, " otscon otslin ConInterOut outlin are%11.4e%11.4e%11.4e%11.4e\n",
01005 rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i],
01006 rfield.outlin[i]+rfield.outlin_noplot[i] );
01007 fprintf( ioQQQ, " line continuum here are %4.4s %4.4s\n",
01008 rfield.chLineLabel[i], rfield.chContLabel[i] );
01009 }
01010
01011 }
01012 }
01013
01014 phisig += rfield.flux[i];
01015
01016
01017 if( phisig > 0. && rfield.SummedDif[i]> 0. )
01018 {
01019 if( fabs(rfield.SummedCon[i]/phisig-1.) > 1e-3 )
01020 {
01021 ++nInsane;
01022
01023
01024 if( nInsane < LIM_INSAME_PRT )
01025 {
01026 fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum %3ld, insane SummedCon at energy %.5e error=%.2e i=%ld\n",
01027 ipPnt, rfield.anu[i], rfield.SummedCon[i]/phisig - 1., i );
01028 fprintf( ioQQQ, " SummedCon, sum are %.4e %.4e\n",
01029 rfield.SummedCon[i], phisig );
01030 fprintf( ioQQQ, " otscon otslin ConInterOut outlin flux are%.4e %.4e %.4e %.4e %.4e\n",
01031 rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i],
01032 rfield.outlin[i]+rfield.outlin_noplot[i], rfield.flux[i] );
01033 fprintf( ioQQQ, " line continuum here are %s %s\n",
01034 rfield.chLineLabel[i], rfield.chContLabel[i]
01035 );
01036 }
01037 }
01038 }
01039 }
01040
01041 if( nInsane > 0 )
01042 {
01043 fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum too much insanity to continue.\n");
01044
01045 TotalInsanity();
01046 }
01047
01048 DEBUG_EXIT( "RT_OTS_ChkSum()" );
01049 return;
01050 }
01051
01052
01053
01054
01055 void RT_OTS_PrtRate(
01056
01057 double weak ,
01058
01059 int chFlag )
01060 {
01061 long int i;
01062
01063 DEBUG_ENTRY( "RT_OTS_PrtRate()" );
01064
01065
01066 ASSERT( chFlag=='l' || chFlag=='c' || chFlag=='b' );
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080 if( chFlag == 'c' || chFlag == 'b' )
01081 {
01082 fprintf( ioQQQ, " DEBUG OTSCON array, anu, otscon, opac, OTS*opac limit:%.2e zone:%.2f IonConv?%c\n",
01083 weak,fnzone ,TorF(conv.lgConvIoniz) );
01084
01085 for( i=0; i < rfield.nupper; i++ )
01086 {
01087 if( rfield.otscon[i]*opac.opacity_abs[i] > weak )
01088 {
01089 fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s \n",
01090 i,
01091 rfield.anu[i],
01092 rfield.otscon[i],
01093 opac.opacity_abs[i],
01094 rfield.otscon[i]*opac.opacity_abs[i],
01095 rfield.chContLabel[i]);
01096
01097 }
01098 }
01099 }
01100
01101
01102
01103
01104 if( chFlag == 'l' || chFlag == 'b' )
01105 {
01106 fprintf( ioQQQ, " DEBUG OTSLIN array, anu, otslin, opac, OTS*opac Lab nLine limit:%.2e zone:%.2f IonConv?%c\n",
01107 weak,fnzone,TorF(conv.lgConvIoniz) );
01108
01109 for( i=0; i < rfield.nupper; i++ )
01110 {
01111 if( rfield.otslin[i]*opac.opacity_abs[i] > weak )
01112 {
01113 fprintf( ioQQQ, " %4ld%12.4e%12.4e%12.4e%12.4e %s %3li\n",
01114 i,
01115 rfield.anu[i],
01116 rfield.otslin[i],
01117 opac.opacity_abs[i],
01118 rfield.otslin[i]*opac.opacity_abs[i],
01119 rfield.chLineLabel[i] ,
01120 rfield.line_count[i] );
01121 }
01122 }
01123 }
01124
01125 DEBUG_EXIT( "RT_OTS_PrtRate()" );
01126 return;
01127 }