00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "optimize.h"
00007 #include "rfield.h"
00008 #include "radius.h"
00009 #include "geometry.h"
00010 #include "iterations.h"
00011 #include "stopcalc.h"
00012 #include "input.h"
00013 #include "parse.h"
00014
00015 void ParseStop(char *chCard)
00016 {
00017 bool lgEOL;
00018
00019 long int i,
00020 j;
00021
00022 double a,
00023 effcol,
00024 tread;
00025
00026 DEBUG_ENTRY( "ParseStop()" );
00027
00028
00029 i = 5;
00030 a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00031
00032
00033 if( lgEOL && !nMatch(" OFF",chCard) )
00034 {
00035 NoNumb(chCard);
00036 }
00037
00038 if( nMatch("TEMP",chCard) )
00039 {
00040
00041 if( lgEOL && nMatch(" OFF",chCard) )
00042 {
00043
00044 StopCalc.tend = -1.f;
00045
00046
00047
00048 }
00049 else
00050 {
00051
00052
00053
00054 if( a <= 10. && !nMatch("LINE",chCard) )
00055 {
00056 tread = pow(10.,a);
00057 }
00058 else
00059 {
00060 tread = a;
00061 }
00062
00063
00064 if( tread < 3. )
00065 {
00066 fprintf( ioQQQ, " Temperatures below 3K not allowed. Reset to 3K. I am doing this myself.\n" );
00067 tread = 3.;
00068 }
00069
00070 if( nMatch("EXCE",chCard) )
00071 {
00072
00073
00074 StopCalc.T2High = (float)tread;
00075 }
00076 else
00077 {
00078
00079 StopCalc.tend = (float)tread;
00080
00081
00082
00083 }
00084 }
00085 }
00086
00087
00088 else if( nMatch("OPTI",chCard) && nMatch("21CM",chCard) )
00089 {
00090
00091 bool lgLOG = true;
00092 if( nMatch("LINE",chCard) )
00093 {
00094
00095 lgLOG = false;
00096 }
00097 i = 4;
00098 j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00099 if( j!=21 )
00100 {
00101 fprintf( ioQQQ, " First number on STOP 21CM OPTICAL DEPTH command must be 21\n" );
00102 puts( "[Stop in ParseStop]" );
00103 cdEXIT(EXIT_FAILURE);
00104 }
00105
00106 a = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00107
00108
00109 if( lgLOG )
00110 {
00111 StopCalc.tauend = (float)pow(10.,a);
00112 }
00113 else
00114 {
00115 StopCalc.tauend = (float)a;
00116 }
00117
00118 StopCalc.lgStop21cm = true;
00119 }
00120
00121 else if( nMatch("OPTI",chCard) )
00122 {
00123
00124 bool lgLOG = true;
00125 if( nMatch("LINE",chCard) )
00126 {
00127
00128 lgLOG = false;
00129 }
00130
00131 if( a > 37. )
00132 {
00133 fprintf( ioQQQ, " optical depth too big\n" );
00134 puts( "[Stop in ParseStop]" );
00135 cdEXIT(EXIT_FAILURE);
00136 }
00137
00138
00139 if( lgLOG )
00140 {
00141 StopCalc.tauend = (float)pow(10.,a);
00142 }
00143 else
00144 {
00145 StopCalc.tauend = (float)a;
00146 }
00147
00148
00149 StopCalc.taunu = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00150
00151 if( lgEOL )
00152 {
00153 if( nMatch("LYMA",chCard) )
00154 {
00155
00156 StopCalc.taunu = 1.;
00157 }
00158 else if( nMatch("BALM",chCard) )
00159 {
00160
00161 StopCalc.taunu = 0.25;
00162 }
00163 else
00164 {
00165 fprintf( ioQQQ, " There must be a second number, the energy in Ryd. Sorry.\n" );
00166 puts( "[Stop in ParseStop]" );
00167 cdEXIT(EXIT_FAILURE);
00168 }
00169 }
00170
00171 else
00172 {
00173
00174 if( StopCalc.taunu < 0. )
00175 {
00176 StopCalc.taunu = (float)pow(10.f,StopCalc.taunu);
00177 }
00178
00179
00180 if( StopCalc.taunu < rfield.emm || StopCalc.taunu > rfield.egamry )
00181 {
00182 fprintf( ioQQQ, " The energy must be in the range %10.2e to %10.2e. It was %10.2e. Sorry.\n",
00183 rfield.emm, rfield.egamry, StopCalc.taunu );
00184 puts( "[Stop in ParseStop]" );
00185 cdEXIT(EXIT_FAILURE);
00186 }
00187 }
00188 }
00189
00190
00191 else if( nMatch(" AV ",chCard) )
00192 {
00193
00194 if( a<=0. )
00195 {
00196 a = pow(10.,a);
00197 }
00198
00199
00200 if( nMatch("EXTE" , chCard ) )
00201 {
00202 StopCalc.AV_extended = (float)a;
00203 }
00204 else
00205 {
00206
00207 StopCalc.AV_point = (float)a;
00208 }
00209 }
00210
00211
00212 else if( nMatch("MOLE",chCard) && nMatch("DEPL",chCard) )
00213 {
00214 if( a <= 0. )
00215 {
00216 StopCalc.StopDepleteFrac = (float)pow(10.,a);
00217 }
00218 else
00219 {
00220 StopCalc.StopDepleteFrac = (float)a;
00221 }
00222 }
00223
00224
00225 else if( nMatch("MASS",chCard) )
00226 {
00227
00228
00229
00230 StopCalc.xMass = (float)a;
00231
00232 if( StopCalc.xMass == 0 )
00233 StopCalc.xMass = SMALLFLOAT;
00234 }
00235
00236
00237
00238 else if( nMatch("THIC",chCard) || nMatch("DEPT",chCard) )
00239 {
00240
00241
00242 if( nMatch("LINE",chCard) )
00243 {
00244 radius.router[0] = a;
00245 }
00246 else
00247 {
00248
00249 if( a > 37. )
00250 {
00251 fprintf( ioQQQ, " thickness too large\n" );
00252 puts( "[Stop in ParseStop]" );
00253 cdEXIT(EXIT_FAILURE);
00254 }
00255 radius.router[0] = pow(10.,a);
00256 }
00257
00258
00259 if( nMatch("PARS",chCard) )
00260 {
00261 radius.router[0] *= PARSEC;
00262 }
00263
00264
00265 for( j=1; j < iterations.iter_malloc; j++ )
00266 {
00267 a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00268 if( lgEOL )
00269 {
00270 radius.router[j] = radius.router[j-1];
00271 }
00272 else
00273 {
00274 if( nMatch("LINE",chCard) )
00275 {
00276 radius.router[j] = a;
00277 }
00278 else
00279 {
00280 if( a > 37. )
00281 {
00282 fprintf( ioQQQ, " thickness too large\n" );
00283 puts( "[Stop in ParseStop]" );
00284 cdEXIT(EXIT_FAILURE);
00285 }
00286 radius.router[j] = pow(10.,a);
00287 }
00288 if( nMatch("PARS",chCard) )
00289 {
00290 radius.router[j] *= PARSEC;
00291 }
00292 }
00293 }
00294
00295
00296 if( optimize.lgVarOn )
00297 {
00298 optimize.nvarxt[optimize.nparm] = 1;
00299 strcpy( optimize.chVarFmt[optimize.nparm], "STOP THICKNESS %f" );
00300
00301
00302 optimize.nvfpnt[optimize.nparm] = input.nRead;
00303
00304
00305 optimize.vparm[0][optimize.nparm] = (float)log10(radius.router[0]);
00306 optimize.vincr[optimize.nparm] = 0.5;
00307
00308 ++optimize.nparm;
00309 }
00310 }
00311
00312
00313 else if( nMatch("ZONE",chCard) )
00314 {
00315
00316
00317
00318 geometry.nend[0] = (long)MAX2(1.,a);
00319 geometry.lgZoneSet = true;
00320
00321
00322 geometry.lgEndDflt = false;
00323
00324 for( j=1; j < iterations.iter_malloc; j++ )
00325 {
00326 geometry.nend[j] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00327
00328
00329 if( lgEOL )
00330 {
00331 geometry.nend[j] = geometry.nend[j-1];
00332 }
00333 else
00334 {
00335
00336
00337 geometry.nend[j] = MAX2( 1 , geometry.nend[j] );
00338 }
00339 }
00340 }
00341
00342
00343 else if( nMatch("EFRA",chCard) )
00344 {
00345 if( a <= 0. )
00346 {
00347 StopCalc.StopElecFrac = (float)pow(10.,a);
00348 }
00349 else
00350 {
00351 StopCalc.StopElecFrac = (float)a;
00352 }
00353 }
00354
00355
00356
00357 else if( nMatch("MFRA",chCard) )
00358 {
00359 if( a <= 0. )
00360 {
00361 StopCalc.StopH2MoleFrac = (float)pow(10.,a);
00362 }
00363 else
00364 {
00365 StopCalc.StopH2MoleFrac = (float)a;
00366 }
00367 }
00368
00369
00370
00371 else if( nMatch("PFRA",chCard) )
00372 {
00373 if( a <= 0. )
00374 {
00375 StopCalc.StopHPlusFrac = (float)pow(10.,a);
00376 }
00377 else
00378 {
00379 StopCalc.StopHPlusFrac = (float)a;
00380 }
00381 }
00382
00383
00384 else if( nMatch("COLU",chCard) )
00385 {
00386
00387 if( nMatch("EFFE",chCard) )
00388 {
00389
00390 effcol = pow(10.,a);
00391 StopCalc.tauend = (float)(effcol*2.14e-22);
00392 StopCalc.taunu = 73.5;
00393
00394 if( optimize.lgVarOn )
00395 {
00396 optimize.nvarxt[optimize.nparm] = 1;
00397 strcpy( optimize.chVarFmt[optimize.nparm], "STOP EFFECTIVE COLUMN DENSITY %f" );
00398
00399 optimize.nvfpnt[optimize.nparm] = input.nRead;
00400
00401 optimize.vparm[0][optimize.nparm] = (float)log10(effcol);
00402 optimize.vincr[optimize.nparm] = 0.5;
00403 ++optimize.nparm;
00404 }
00405 }
00406
00407 else if( nMatch("IONI",chCard) )
00408 {
00409
00410 if( a > 37. )
00411 {
00412 fprintf( ioQQQ, " column too big\n" );
00413 puts( "[Stop in ParseStop]" );
00414 cdEXIT(EXIT_FAILURE);
00415 }
00416
00417 StopCalc.colpls = (float)pow(10.,a);
00418
00419
00420 if( optimize.lgVarOn )
00421 {
00422 optimize.nvarxt[optimize.nparm] = 1;
00423 strcpy( optimize.chVarFmt[optimize.nparm], "STOP IONIZED COLUMN DENSITY %f" );
00424
00425 optimize.nvfpnt[optimize.nparm] = input.nRead;
00426
00427 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.colpls);
00428 optimize.vincr[optimize.nparm] = 0.5;
00429 ++optimize.nparm;
00430 }
00431 }
00432
00433
00434 else if( nMatch("NEUT",chCard) )
00435 {
00436 StopCalc.colnut = (float)pow(10.,a);
00437
00438
00439 if( optimize.lgVarOn )
00440 {
00441 optimize.nvarxt[optimize.nparm] = 1;
00442 strcpy( optimize.chVarFmt[optimize.nparm], "STOP NEUTRAL COLUMN DENSITY %f");
00443
00444 optimize.nvfpnt[optimize.nparm] = input.nRead;
00445
00446 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.colnut);
00447 optimize.vincr[optimize.nparm] = 0.5;
00448 ++optimize.nparm;
00449 }
00450 }
00451
00452
00453
00454
00455 else if( nMatch(" H2 ",chCard) )
00456 {
00457
00458
00459
00460 i = 5;
00461 j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00462 if( j != 2 )
00463 {
00464 fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
00465 fprintf( ioQQQ, " The first number I encounter should be the 2 in H2.\n Sorry.\n" );
00466 puts( "[Stop in ParseStop]" );
00467 cdEXIT(EXIT_FAILURE);
00468 }
00469 a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00470 StopCalc.col_h2 = (float)pow(10.,a);
00471
00472
00473 if( optimize.lgVarOn )
00474 {
00475 optimize.nvarxt[optimize.nparm] = 1;
00476 strcpy( optimize.chVarFmt[optimize.nparm], "STOP H2 COLUMN DENSITY %f");
00477
00478 optimize.nvfpnt[optimize.nparm] = input.nRead;
00479
00480 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_h2);
00481 optimize.vincr[optimize.nparm] = 0.5;
00482 ++optimize.nparm;
00483 }
00484 }
00485
00486 else if( nMatch("ATOM",chCard) )
00487 {
00488 StopCalc.col_h2_nut = (float)pow(10.,a);
00489
00490 if( optimize.lgVarOn )
00491 {
00492 optimize.nvarxt[optimize.nparm] = 1;
00493 strcpy( optimize.chVarFmt[optimize.nparm], "STOP ATOMIC COLUMN DENSITY %f");
00494
00495 optimize.nvfpnt[optimize.nparm] = input.nRead;
00496
00497 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_h2_nut);
00498 optimize.vincr[optimize.nparm] = 0.5;
00499 ++optimize.nparm;
00500 }
00501 }
00502
00503 else if( nMatch("H/TS",chCard) )
00504 {
00505
00506 StopCalc.col_H0_ov_Tspin = (float)pow(10.,a);
00507
00508 if( optimize.lgVarOn )
00509 {
00510 optimize.nvarxt[optimize.nparm] = 1;
00511 strcpy( optimize.chVarFmt[optimize.nparm], "STOP H/TSPIN COLUMN DENSITY %f");
00512
00513 optimize.nvfpnt[optimize.nparm] = input.nRead;
00514
00515 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_H0_ov_Tspin);
00516 optimize.vincr[optimize.nparm] = 0.5;
00517 ++optimize.nparm;
00518 }
00519 }
00520
00521 else if( nMatch(" CO ",chCard) )
00522 {
00523
00524
00525 StopCalc.col_monoxco = (float)pow(10.,a);
00526
00527 if( optimize.lgVarOn )
00528 {
00529 optimize.nvarxt[optimize.nparm] = 1;
00530 strcpy( optimize.chVarFmt[optimize.nparm], "STOP CO COLUMN DENSITY %f");
00531
00532 optimize.nvfpnt[optimize.nparm] = input.nRead;
00533
00534 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.col_monoxco);
00535 optimize.vincr[optimize.nparm] = 0.5;
00536 ++optimize.nparm;
00537 }
00538 }
00539
00540
00541 else
00542 {
00543
00544 if( a > 37. )
00545 {
00546 fprintf( ioQQQ, " column too big\n" );
00547 puts( "[Stop in ParseStop]" );
00548 cdEXIT(EXIT_FAILURE);
00549 }
00550
00551 StopCalc.HColStop = (float)pow(10.,a);
00552
00553
00554 if( optimize.lgVarOn )
00555 {
00556 optimize.nvarxt[optimize.nparm] = 1;
00557 strcpy( optimize.chVarFmt[optimize.nparm], "STOP COLUMN DENSITY %f" );
00558
00559 optimize.nvfpnt[optimize.nparm] = input.nRead;
00560
00561 optimize.vparm[0][optimize.nparm] = (float)log10(StopCalc.HColStop);
00562 optimize.vincr[optimize.nparm] = 0.5;
00563 ++optimize.nparm;
00564 }
00565 }
00566 }
00567
00568
00569 else if( nMatch("EDEN",chCard) )
00570 {
00571
00572
00573 if( nMatch("LINE",chCard) )
00574 {
00575 StopCalc.StopElecDensity = (float)a;
00576 }
00577 else
00578 {
00579 StopCalc.StopElecDensity = (float)pow(10.,a);
00580 }
00581 }
00582
00583
00584
00585 else if( nMatch("LINE",chCard) )
00586 {
00587 char chLabel[5];
00588
00589
00590
00591
00592
00593 GetQuote( chLabel , chCard , true );
00594
00595
00596 cap4( StopCalc.chStopLabel[StopCalc.nstpl], chLabel);
00597
00598
00599
00600 i = 5;
00601
00602 StopCalc.StopLineWl[StopCalc.nstpl] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
00603
00604
00605 if( input.chCARDCAPS[i-1] == 'M' )
00606 {
00607
00608 StopCalc.StopLineWl[StopCalc.nstpl] *= 1e4;
00609 }
00610 else if( input.chCARDCAPS[i-1] == 'C' )
00611 {
00612
00613 StopCalc.StopLineWl[StopCalc.nstpl] *= 1e8;
00614 }
00615
00616
00617 a = log10( StopCalc.StopLineWl[StopCalc.nstpl] );
00618 if( StopCalc.StopLineWl[StopCalc.nstpl]>=1. )
00619 {
00620 a -= floor(a);
00621 }
00622 else
00623 {
00624 a += floor(a);
00625 }
00626 a = pow(10.,a);
00627
00628 StopCalc.ErrorLineWl[StopCalc.nstpl] = (float)(0.00105/a);
00629
00630
00631 StopCalc.stpint[StopCalc.nstpl] =
00632 (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
00633
00634 if( lgEOL )
00635 {
00636 fprintf( ioQQQ, " There MUST be a relative intensity entered. Sorry.\n" );
00637 puts( "[Stop in ParseStop]" );
00638 cdEXIT(EXIT_FAILURE);
00639 }
00640
00641
00642
00643
00644 j = i;
00645 a = FFmtRead(chCard,&j,INPUT_LINE_LENGTH, &lgEOL);
00646
00647 if( lgEOL )
00648 {
00649
00650 StopCalc.StopLineWl2[StopCalc.nstpl] = 1;
00651 }
00652 else
00653 {
00654
00655
00656 GetQuote( chLabel , chCard , true );
00657
00658
00659 cap4( StopCalc.chStopLabel2[StopCalc.nstpl], chLabel);
00660
00661
00662
00663 StopCalc.StopLineWl2[StopCalc.nstpl] =
00664 (float)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);
00665
00666
00667 if( input.chCARDCAPS[i-1] == 'M' )
00668 {
00669
00670 StopCalc.StopLineWl2[StopCalc.nstpl] *= 1e4;
00671 }
00672 else if( input.chCARDCAPS[i-1] == 'C' )
00673 {
00674
00675 StopCalc.StopLineWl2[StopCalc.nstpl] *= 1e8;
00676 }
00677
00678
00679 if( StopCalc.StopLineWl2[StopCalc.nstpl] > 0. )
00680 {
00681 a = log10( StopCalc.StopLineWl2[StopCalc.nstpl] );
00682 if( StopCalc.StopLineWl2[StopCalc.nstpl]>=1. )
00683 {
00684 a -= floor(a);
00685 }
00686 else
00687 {
00688 a += floor(a);
00689 }
00690 a = pow(10.,a);
00691
00692 StopCalc.ErrorLineWl2[StopCalc.nstpl] = (float)(0.00105/a);
00693 }
00694 else
00695 {
00696 StopCalc.ErrorLineWl2[StopCalc.nstpl] = 1e-4f;
00697 }
00698 }
00699
00700
00701 StopCalc.nstpl = MIN2(StopCalc.nstpl+1,MXSTPL-1);
00702 }
00703
00704 else if( nMatch("NTOT" , chCard ) && nMatch("ALIO" , chCard ) )
00705 {
00706
00707
00708
00709 StopCalc.nTotalIonizStop = (long)a;
00710 }
00711
00712
00713 else
00714 {
00715 fprintf( ioQQQ, " keyword error for STOP line, line image follows;\n" );
00716 fprintf( ioQQQ, "==%s==\n" , chCard);
00717 puts( "[Stop in ParseStop]" );
00718 cdEXIT(EXIT_FAILURE);
00719 }
00720
00721 DEBUG_EXIT( "ParseStop()" );
00722 return;
00723 }
00724