00001
00002
00003
00004 #include "cddefines.h"
00005 #include "elementnames.h"
00006 #include "optimize.h"
00007 #include "hydrogenic.h"
00008 #include "input.h"
00009 #include "iso.h"
00010 #include "parser.h"
00011 #include "phycon.h"
00012 #include "rfield.h"
00013 #include "taulines.h"
00014 #include "thirdparty.h"
00015
00016
00017 void ParseAtomISO(long ipISO, Parser &p )
00018 {
00019 long int numLevels;
00020
00021 DEBUG_ENTRY( "ParseAtomISO()" );
00022
00023
00024
00025 long int nelem = p.GetElem( );
00026
00027
00028 if( ipISO==ipHE_LIKE && nelem==ipHYDROGEN )
00029 {
00030 fprintf(ioQQQ," Sorry, H-like He does not exist.\n");
00031 cdEXIT(EXIT_FAILURE);
00032 }
00033
00034
00035 if( p.nMatch("COLL") && !p.nMatch("LEVE" ) )
00036 {
00037
00038
00039 if( p.nMatch("EXCI") )
00040 {
00041
00042 iso_ctrl.lgColl_excite[ipISO] = false;
00043 phycon.lgPhysOK = false;
00044 }
00045 else if( p.nMatch("IONI") )
00046 {
00047
00048 iso_ctrl.lgColl_ionize[ipISO] = false;
00049 phycon.lgPhysOK = false;
00050 }
00051
00052 else if( p.nMatch("2S2P") || ( p.nMatch("2P2S") && ipISO == ipH_LIKE ) )
00053 {
00054
00055
00056 fprintf(ioQQQ,"This command changed to ATOM H-LIKE COLLISIONS L-MIXING\n");
00057 fprintf(ioQQQ,"I will parse it for now, but may not in the future.\n");
00058
00059 iso_ctrl.lgColl_l_mixing[ipISO] = false;
00060 phycon.lgPhysOK = false;
00061 }
00062
00063 else if( p.nMatch("L-MI") )
00064 {
00065 if( ipISO == ipH_LIKE )
00066 {
00067
00068
00069 iso_ctrl.lgColl_l_mixing[ipISO] = false;
00070 phycon.lgPhysOK = false;
00071 }
00072 else if( p.nMatch("THER") )
00073 {
00074
00075
00076 if( p.nMatch("NO T") )
00077 {
00078
00079
00080
00081 iso_ctrl.lgCS_therm_ave[ipISO] = false;
00082 }
00083 else
00084 {
00085 iso_ctrl.lgCS_therm_ave[ipISO] = true;
00086 }
00087 }
00088 else if( p.nMatch("PENG") )
00089 {
00090 iso_ctrl.lgCS_Vrinceanu[ipISO] = false;
00091 }
00092 else if( p.nMatch(" OFF" ) )
00093 {
00094
00095
00096 iso_ctrl.lgColl_l_mixing[ipISO] = false;
00097 phycon.lgPhysOK = false;
00098 iso_ctrl.lgCS_Vrinceanu[ipISO] = false;
00099 iso_ctrl.lgCS_therm_ave[ipISO] = false;
00100 }
00101 else
00102 {
00103 fprintf( ioQQQ, " needs parameter\n" );
00104 cdEXIT(EXIT_FAILURE);
00105 }
00106 }
00107 else if( p.nMatch(" OFF" ) )
00108 {
00109
00110 iso_ctrl.lgColl_excite[ipISO] = false;
00111 iso_ctrl.lgColl_ionize[ipISO] = false;
00112 iso_ctrl.lgColl_l_mixing[ipISO] = false;
00113 phycon.lgPhysOK = false;
00114 }
00115 else
00116 {
00117 fprintf( ioQQQ, " needs parameter\n" );
00118 cdEXIT(EXIT_FAILURE);
00119 }
00120 }
00121
00122 else if( p.nMatch("CONT") && p.nMatch("LOWE") )
00123 {
00124
00125 if( p.nMatch("OFF") )
00126 iso_ctrl.lgContinuumLoweringEnabled[ipISO] = false;
00127 else
00128 iso_ctrl.lgContinuumLoweringEnabled[ipISO] = true;
00129 }
00130
00131 else if( p.nMatch("DAMP") )
00132 {
00133 if( ipISO == ipHE_LIKE )
00134 {
00135 fprintf(ioQQQ," Sorry, the DAMPING option is not implemented for the he-like sequence.\n");
00136 cdEXIT(EXIT_FAILURE);
00137 }
00138
00139
00140 hydro.DampOnFac = 0.;
00141 }
00142
00143 else if( p.nMatch("DIEL") )
00144 {
00145 if( ipISO == ipH_LIKE )
00146 {
00147 fprintf(ioQQQ," Sorry, but dielectronic recombination onto the h-like sequence is not possible.\n");
00148 cdEXIT(EXIT_FAILURE);
00149 }
00150
00151
00152 if( p.nMatch(" OFF") )
00153 {
00154 iso_ctrl.lgDielRecom[ipISO] = false;
00155 }
00156 else
00157 iso_ctrl.lgDielRecom[ipISO] = true;
00158 }
00159
00160 else if( p.nMatch("LEVE") )
00161 {
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 if( p.nMatch("LTE") )
00175 {
00176
00177 iso_ctrl.lgLTE_levels[ipISO] = true;
00178 }
00179 else if( p.nMatch("PRIN") )
00180 {
00181
00182 iso_ctrl.lgPrintNumberOfLevels = true;
00183 }
00184 else if( !lgHydroMalloc )
00185 {
00186 numLevels = (long int)p.FFmtRead();
00187
00188 if( !p.lgEOL() )
00189 {
00190 if( ipISO == ipH_LIKE && numLevels > NHYDRO_MAX_LEVEL-2 )
00191 {
00192 fprintf( ioQQQ, " Not possible to set nhlvl to >NHYDRO_MAX_LEVEL-2= %i\n",
00193 NHYDRO_MAX_LEVEL-2 );
00194 fprintf( ioQQQ, " change NHYDRO_MAX_LEVEL\n");
00195 cdEXIT(EXIT_FAILURE);
00196 }
00197
00198
00199 if( !p.nMatch("COLL") && ipISO == ipH_LIKE &&
00200 ( 2. / POW3((double)numLevels) < rfield.emm ) )
00201 {
00202 fprintf( ioQQQ, " Not possible to set iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_max to such a high value, since "
00203 "alpha transition not within energy bounds of code\n");
00204
00205 fprintf( ioQQQ, " lowest energy is %e and corresponding highest level is %li\n" ,
00206 rfield.emm, (long)pow(2./rfield.emm, 0.3333) );
00207 cdEXIT(EXIT_FAILURE);
00208 }
00209 }
00210
00211 if( p.lgEOL() )
00212 {
00213 int LevelsResolved=-1 , LevelsCollapsed=10;
00214
00215 if( p.nMatch("LARG") )
00216 {
00217
00218 LevelsResolved = RREC_MAXN;
00219 }
00220
00221
00222 else if( p.nMatch("SMAL") || p.nMatch("COMP") )
00223 {
00224 if( ipISO == ipH_LIKE )
00225 LevelsResolved = 5;
00226 else if( ipISO == ipHE_LIKE )
00227 LevelsResolved = 3;
00228 }
00229 else
00230
00231 p.NoNumb("levels");
00232
00233 if( nelem<0 )
00234 {
00235
00236 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00237 {
00238 iso_sp[ipISO][nelem].nCollapsed_max =
00239 MIN2( iso_sp[ipISO][nelem].nCollapsed_max , LevelsCollapsed );
00240 iso_sp[ipISO][nelem].n_HighestResolved_max =
00241 MIN2( iso_sp[ipISO][nelem].n_HighestResolved_max , LevelsResolved );
00242 iso_update_num_levels( ipISO, nelem );
00243 }
00244 }
00245 else
00246 {
00247 iso_sp[ipISO][nelem].nCollapsed_max = LevelsCollapsed;
00248 iso_sp[ipISO][nelem].n_HighestResolved_max = LevelsResolved;
00249 iso_update_num_levels( ipISO, nelem );
00250 }
00251 }
00252
00253 else if( p.nMatch("COLLAP") )
00254 {
00255
00256 if( numLevels < 1 )
00257 {
00258 fprintf( ioQQQ, "There must be at least one collapsed level.\n");
00259 cdEXIT(EXIT_FAILURE);
00260 }
00261
00262 if( nelem<0 )
00263 {
00264
00265 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00266 {
00267 iso_sp[ipISO][nelem].nCollapsed_max = numLevels;
00268 iso_update_num_levels( ipISO, nelem );
00269 }
00270 }
00271 else
00272 {
00273 iso_sp[ipISO][nelem].nCollapsed_max = numLevels;
00274 iso_update_num_levels( ipISO, nelem );
00275 }
00276 }
00277 else if( p.nMatch("RESOLV") )
00278 {
00279
00280 if( ( numLevels < 3 ) && !p.nMatch("COLL") )
00281 {
00282 fprintf( ioQQQ, " cannot have fewer than 3 resolved levels, the requested number was %li\n" ,
00283 numLevels );
00284 fprintf( ioQQQ, " Sorry.\n" );
00285 cdEXIT(EXIT_FAILURE);
00286 }
00287
00288 if( nelem<0 )
00289 {
00290
00291 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00292 {
00293 iso_sp[ipISO][nelem].n_HighestResolved_max = numLevels;
00294 iso_update_num_levels( ipISO, nelem );
00295 }
00296 }
00297 else
00298 {
00299 iso_sp[ipISO][nelem].n_HighestResolved_max = numLevels;
00300 iso_update_num_levels( ipISO, nelem );
00301 }
00302 }
00303 else
00304 {
00305 fprintf(ioQQQ, "I did not recognize a keyword on this atom xx-like levels command."
00306 " Should be COLLAPSED or RESOLVED.\n Sorry.\n\n");
00307 cdEXIT(EXIT_FAILURE);
00308 }
00309 }
00310 }
00311
00312 else if( p.nMatch("ERRO") && p.nMatch("GENE" ) )
00313 {
00314
00315
00316 iso_ctrl.lgRandErrGen[ipISO] = true;
00317 iso_ctrl.modelRank[ipISO] = (int)p.FFmtRead();
00318
00319 iso_ctrl.modelRank[ipISO] = MAX2( 0, iso_ctrl.modelRank[ipISO] );
00320 if( p.lgEOL() )
00321
00322
00323 iso_ctrl.modelRank[ipISO] = abs((int)time(NULL));
00324
00325
00326
00327
00328 init_genrand( (unsigned)iso_ctrl.modelRank[ipISO] + 2);
00329
00330 if( p.nMatch("PESS") )
00331 iso_ctrl.lgPessimisticErrors = true;
00332 else
00333 iso_ctrl.lgPessimisticErrors = false;
00334 }
00335
00336 else if( p.nMatch(" FSM") )
00337 {
00338 if( ipISO == ipH_LIKE )
00339 {
00340 fprintf(ioQQQ," Sorry, but fine-structure mixing can only be implemented for the He-like sequence.\n");
00341 cdEXIT(EXIT_FAILURE);
00342 }
00343
00344
00345
00346 if( p.nMatch(" OFF") )
00347 iso_ctrl.lgFSM[ipISO] = false;
00348 else
00349 iso_ctrl.lgFSM[ipISO] = true;
00350 }
00351
00352 else if( p.nMatch("GBAR") )
00353 {
00354 if( ipISO == ipH_LIKE )
00355 {
00356 fprintf(ioQQQ," Sorry, the GBAR option is only implemented for the He-like sequence.\n");
00357 cdEXIT(EXIT_FAILURE);
00358 }
00359
00360
00361
00362 iso_ctrl.lgCS_Vriens[ipISO] = false;
00363 iso_ctrl.lgCS_None[ipISO] = false;
00364 iso_ctrl.nCS_new[ipISO] = false;
00365
00366
00367 if( p.nMatch("VRIE") )
00368 {
00369
00370 iso_ctrl.lgCS_Vriens[ipISO] = true;
00371 }
00372 else if( p.nMatch(" NEW") )
00373 {
00374
00375 iso_ctrl.nCS_new[ipISO] = (int)p.FFmtRead();
00376
00377
00378
00379 if( p.lgEOL() )
00380 iso_ctrl.nCS_new[ipISO] = 1;
00381
00382 ASSERT( iso_ctrl.nCS_new[ipISO] );
00383 }
00384 else if( p.nMatch(" OFF") )
00385 {
00386
00387 iso_ctrl.lgCS_None[ipISO] = true;
00388 }
00389 else
00390 {
00391 fprintf( ioQQQ, " needs parameter\n" );
00392 cdEXIT(EXIT_FAILURE);
00393 }
00394 }
00395
00396
00397 else if( p.nMatch("LYMA") )
00398 {
00399 if( ipISO == ipH_LIKE && p.nMatch("PUMP") )
00400 {
00401
00402 if( p.nMatch(" OFF") )
00403 {
00404
00405 hydro.lgLymanPumping = false;
00406 }
00407 else if( p.nMatch("SCALE") )
00408 {
00409
00410
00411
00412 hydro.xLymanPumpingScaleFactor =
00413 (realnum)p.FFmtRead();
00414
00415
00416
00417 if( hydro.xLymanPumpingScaleFactor <= 0. ||
00418 p.nMatch(" LOG") )
00419 {
00420 hydro.xLymanPumpingScaleFactor =
00421 (realnum)pow((realnum)10.f , hydro.xLymanPumpingScaleFactor );
00422 }
00423
00424
00425 if( optimize.lgVarOn )
00426 {
00427 optimize.nvarxt[optimize.nparm] = 1;
00428 strcpy( optimize.chVarFmt[optimize.nparm], "ATOM H-LIKE LYMAN PUMPING SCALE %f LOG" );
00429
00430
00431 optimize.nvfpnt[optimize.nparm] = input.nRead;
00432
00433
00434 optimize.vincr[optimize.nparm] = 0.1f;
00435 optimize.vparm[0][optimize.nparm] = (realnum)log10(hydro.xLymanPumpingScaleFactor);
00436 ++optimize.nparm;
00437 }
00438 }
00439 else
00440 {
00441 fprintf(ioQQQ," Sorry, I didn\'t recognize an option on this ATOM H-LIKE LYMAN PUMP command.\n");
00442 fprintf(ioQQQ," The options are \" OFF\", and \"SCALE\".\n");
00443 cdEXIT(EXIT_FAILURE);
00444 }
00445 }
00446 else if( p.nMatch("EXTRA") )
00447 {
00448
00449 iso_ctrl.nLyman[ipISO] = (long int)p.FFmtRead();
00450 if( p.lgEOL() )
00451 p.NoNumb("'extra' Lyman lines");
00452 }
00453 else
00454 {
00455 fprintf(ioQQQ," Sorry, I didn\'t recognize an option on this ATOM xx-LIKE LYMAN command.\n");
00456 fprintf(ioQQQ," The options are \"PUMP\", and \"EXTRA\".\n");
00457 cdEXIT(EXIT_FAILURE);
00458 }
00459 }
00460
00461
00462 else if( p.nMatch("RECO") &&
00463 p.nMatch(" NO ") && p.nMatch("INTE") )
00464 {
00465
00466
00467
00468 iso_ctrl.lgNoRecombInterp[ipISO] = true;
00469 }
00470
00471 else if( p.nMatch("REDI") )
00472 {
00473 int ipRedis=0;
00474
00475
00476
00477
00478
00479 if( p.nMatch(" PRD") )
00480 {
00481 ipRedis = ipPRD;
00482 }
00483
00484 else if( p.nMatch(" CRD") )
00485 {
00486 ipRedis = ipCRD;
00487 }
00488
00489 else if( p.nMatch("CRDW") )
00490 {
00491 ipRedis = ipCRDW;
00492 }
00493
00494
00495 else if( !p.nMatch("SHOW") )
00496 {
00497 fprintf(ioQQQ," There should have been a second keyword on this command.\n");
00498 fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n");
00499 cdEXIT(EXIT_FAILURE);
00500 }
00501
00502
00503 if( p.nMatch("ALPH") )
00504 {
00505 iso_ctrl.ipLyaRedist[ipISO] = ipRedis;
00506 }
00507
00508 else if( p.nMatch("RESO") )
00509 {
00510 iso_ctrl.ipResoRedist[ipISO] = ipRedis;
00511 }
00512
00513 else if( p.nMatch("SUBO") )
00514 {
00515 iso_ctrl.ipSubRedist[ipISO] = ipRedis;
00516 }
00517
00518 else if( p.nMatch("SHOW") )
00519 {
00520 fprintf(ioQQQ," Ly a is ");
00521 if( iso_ctrl.ipLyaRedist[ipISO] ==ipCRDW )
00522 {
00523 fprintf(ioQQQ,"complete redistribution with wings\n");
00524 }
00525 else if( iso_ctrl.ipLyaRedist[ipISO] ==ipCRD )
00526 {
00527 fprintf(ioQQQ,"complete redistribution with core only.\n");
00528 }
00529 else if( iso_ctrl.ipLyaRedist[ipISO] ==ipPRD )
00530 {
00531 fprintf(ioQQQ,"partial redistribution.\n");
00532 }
00533 else if( iso_ctrl.ipLyaRedist[ipISO] ==ipLY_A )
00534 {
00535 fprintf(ioQQQ,"special Lya.\n");
00536 }
00537 else
00538 {
00539 fprintf(ioQQQ," PROBLEM Impossible value for iso_ctrl.ipLyaRedist.\n");
00540 TotalInsanity();
00541 }
00542
00543 fprintf(ioQQQ," Other %s resonance lines are ",
00544 elementnames.chElementSym[ipISO] );
00545
00546 if( iso_ctrl.ipResoRedist[ipISO] ==ipCRDW )
00547 {
00548 fprintf(ioQQQ,"complete redistribution with wings\n");
00549 }
00550 else if( iso_ctrl.ipResoRedist[ipISO] ==ipCRD )
00551 {
00552 fprintf(ioQQQ,"complete redistribution with core only.\n");
00553 }
00554 else if( iso_ctrl.ipResoRedist[ipISO] ==ipPRD )
00555 {
00556 fprintf(ioQQQ,"partial redistribution.\n");
00557 }
00558 else
00559 {
00560 fprintf(ioQQQ," PROBLEM Impossible value for iso_ctrl.ipResoRedist.\n");
00561 TotalInsanity();
00562 }
00563
00564 fprintf(ioQQQ," %s subordinate lines are ",
00565 elementnames.chElementSym[ipISO] );
00566
00567 if( iso_ctrl.ipSubRedist[ipISO] ==ipCRDW )
00568 {
00569 fprintf(ioQQQ,"complete redistribution with wings\n");
00570 }
00571 else if( iso_ctrl.ipSubRedist[ipISO] ==ipCRD )
00572 {
00573 fprintf(ioQQQ,"complete redistribution with core only.\n");
00574 }
00575 else if( iso_ctrl.ipSubRedist[ipISO] ==ipPRD )
00576 {
00577 fprintf(ioQQQ,"partial redistribution.\n");
00578 }
00579 else
00580 {
00581 fprintf(ioQQQ," PROBLEM Impossible value for iso_ctrl.ipSubRedist.\n");
00582 TotalInsanity();
00583 }
00584 }
00585 else
00586 {
00587 fprintf(ioQQQ," here should have been another keyword on this command.\n");
00588 fprintf(ioQQQ," Options are ALPHA, RESONANCE, SUBORDINATE. Sorry.\n");
00589 cdEXIT(EXIT_FAILURE);
00590 }
00591 }
00592
00593 else if( p.nMatch("TOPO") )
00594 {
00595 if( p.nMatch(" OFF") )
00596 {
00597 iso_ctrl.lgTopoff[ipISO] = false;
00598 fprintf( ioQQQ, "ISO %li TOPOFF is OFF\n", ipISO );
00599 }
00600 else
00601 iso_ctrl.lgTopoff[ipISO] = true;
00602 }
00603
00604
00605 else
00606 {
00607 fprintf( ioQQQ, " There should have been a keyword on this ATOM H-LIKE or HE-LIKE command.\n Sorry.\n" );
00608 cdEXIT(EXIT_FAILURE);
00609 }
00610 return;
00611 }