00001
00002
00003 #include "cddefines.h"
00004 #include "lines_service.h"
00005 #include "physconst.h"
00006 #include "taulines.h"
00007 #include "trace.h"
00008 #include "string.h"
00009 #include "thirdparty.h"
00010 #include <fstream>
00011
00012 static const bool DEBUGSTATE = false;
00013
00014 void atmdat_CHIANTI_readin( long intNS, char *chPrefix )
00015 {
00016 DEBUG_ENTRY( "atmdat_CHIANTI_readin()" );
00017
00018 int intsplinepts,intTranType,intxs;
00019 long int nMolLevs;
00020 FILE *ioElecCollData=NULL, *ioProtCollData=NULL;
00021 realnum fstatwt,fenergyK,fenergyWN,fWLAng,fenergy,feinsteina;
00022 double fScalingParam,fEnergyDiff,fGF,*xs,*spl,*spl2;
00023
00024 char chLine[FILENAME_PATH_LENGTH_2] ,
00025 chEnFilename[FILENAME_PATH_LENGTH_2],
00026 chTraFilename[FILENAME_PATH_LENGTH_2],
00027 chEleColFilename[FILENAME_PATH_LENGTH_2],
00028 chProColFilename[FILENAME_PATH_LENGTH_2];
00029
00030 realnum *dBaseStatesEnergy;
00031 long int *intNewIndex=NULL,*intOldIndex=NULL;
00032 int interror;
00033 bool lgProtonData=false,lgEneLevOrder;
00034
00035
00036 static const int MAX_NUM_LEVELS = 999;
00037
00038 Species[intNS].lgMolecular = false;
00039
00040 strcpy( chEnFilename , chPrefix );
00041 strcpy( chTraFilename , chPrefix );
00042 strcpy( chEleColFilename , chPrefix );
00043 strcpy( chProColFilename , chPrefix );
00044
00045
00046
00047 strcat( chEnFilename , ".elvlc");
00048 uncaps( chEnFilename );
00049 if(DEBUGSTATE)
00050 printf("The data file is %s \n",chEnFilename);
00051
00052
00053 if( trace.lgTrace )
00054 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEnFilename);
00055
00056 fstream elvlcstream;
00057 open_data(elvlcstream, chEnFilename,mode_r);
00058
00059
00060 strcat( chTraFilename , ".wgfa");
00061 uncaps( chTraFilename );
00062 if(DEBUGSTATE)
00063 printf("The data file is %s \n",chTraFilename);
00064
00065
00066 if( trace.lgTrace )
00067 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chTraFilename);
00068
00069 fstream wgfastream;
00070 open_data(wgfastream, chTraFilename,mode_r);
00071
00072
00073 strcat( chEleColFilename , ".splups");
00074 uncaps( chEleColFilename );
00075 if(DEBUGSTATE)
00076 printf("The data file is %s \n",chEleColFilename);
00077
00078 if( trace.lgTrace )
00079 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chEleColFilename);
00080
00081 ioElecCollData = open_data( chEleColFilename, "r" );
00082
00083
00084 strcat( chProColFilename , ".psplups");
00085 uncaps( chProColFilename );
00086 if(DEBUGSTATE)
00087 printf("The data file is %s \n",chProColFilename);
00088
00089 if( trace.lgTrace )
00090 fprintf( ioQQQ," atmdat_CHIANTI_readin opening %s:",chProColFilename);
00091
00092
00093 if( ( ioProtCollData = fopen( chProColFilename , "r" ) ) != NULL )
00094 {
00095 lgProtonData = true;
00096
00097
00098 }
00099 else
00100 {
00101 lgProtonData = false;
00102 }
00103
00104
00105
00106 const int eof_col = 5;
00107 nMolLevs = -1;
00108 lgEneLevOrder = true;
00109 if (elvlcstream.is_open())
00110 {
00111 int nj = 0;
00112 char otemp[eof_col];
00113
00114 while(nj != -1)
00115 {
00116 elvlcstream.get(otemp,eof_col);
00117 elvlcstream.ignore(INT_MAX,'\n');
00118 nj = atoi(otemp);
00119 nMolLevs++;
00120 }
00121 elvlcstream.seekg(0,ios::beg);
00122 }
00123
00124 long HighestIndexInFile = nMolLevs;
00125 nMolLevs = MIN2( nMolLevs, MAX_NUM_LEVELS );
00126
00127 if( nMolLevs <= 0 )
00128 {
00129 fprintf( ioQQQ, "The number of energy levels is non-positive in datafile %s.\n", chEnFilename );
00130 fprintf( ioQQQ, "The file must be corrupted.\n" );
00131 cdEXIT( EXIT_FAILURE );
00132 }
00133
00134 Species[intNS].numLevels_max = nMolLevs;
00135 Species[intNS].numLevels_local = Species[intNS].numLevels_max;
00136
00137
00138 dBaseStatesEnergy = (realnum*)MALLOC((unsigned long)(nMolLevs)*sizeof(realnum));
00139
00140 dBaseStates[intNS] = (quantumState *)MALLOC((size_t)(nMolLevs)*sizeof(quantumState));
00141
00142 dBaseTrans[intNS] = (transition **)MALLOC(sizeof(transition*)*(unsigned long)nMolLevs);
00143 dBaseTrans[intNS][0] = NULL;
00144 for( long ipHi = 1; ipHi < nMolLevs; ipHi++)
00145 {
00146 dBaseTrans[intNS][ipHi] = (transition *)MALLOC(sizeof(transition)*(unsigned long)ipHi);
00147 for( long ipLo = 0; ipLo < ipHi; ipLo++)
00148 {
00149 dBaseTrans[intNS][ipHi][ipLo].Junk();
00150 dBaseTrans[intNS][ipHi][ipLo].Lo = &dBaseStates[intNS][ipLo];
00151 dBaseTrans[intNS][ipHi][ipLo].Hi = &dBaseStates[intNS][ipHi];
00152 }
00153 }
00154
00155 for( long ipLev = 0; ipLev<NUM_COLLIDERS; ipLev++ )
00156 {
00157 CollRatesArray[intNS][ipLev] = NULL;
00158 }
00159
00160 CollRatesArray[intNS][0] = (double**)MALLOC((nMolLevs)*sizeof(double*));
00161 for( long ipHi = 0; ipHi<nMolLevs; ipHi++ )
00162 {
00163 CollRatesArray[intNS][0][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
00164 for( long ipLo = 0; ipLo<nMolLevs; ipLo++)
00165 {
00166 CollRatesArray[intNS][0][ipHi][ipLo] = 0.0;
00167 }
00168 }
00169
00170 if(lgProtonData)
00171 {
00172 CollRatesArray[intNS][1] = (double**)MALLOC((nMolLevs)*sizeof(double*));
00173 for( long ipHi = 0; ipHi<nMolLevs; ipHi++ )
00174 {
00175 CollRatesArray[intNS][1][ipHi] = (double*)MALLOC((unsigned long)(nMolLevs)*sizeof(double));
00176 for( long ipLo = 0; ipLo<nMolLevs; ipLo++)
00177 {
00178 CollRatesArray[intNS][1][ipHi][ipLo] = 0.0;
00179 }
00180 }
00181 }
00182
00183
00184
00185
00186 const int lvl_nrg_col=16;
00187
00188 const int lvl_skipto_nrg = 40;
00189
00190 const int lvl_skip_ryd = 15;
00191
00192
00193 for( long ipLev=0; ipLev<nMolLevs; ipLev++ )
00194 {
00195 if(elvlcstream.is_open())
00196 {
00197 char thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
00198 elvlcstream.seekg(lvl_skipto_nrg,ios::cur);
00199 elvlcstream.get(thtemp,lvl_nrg_col);
00200 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
00201 fenergy = (realnum) atof(thtemp);
00202
00203
00204
00205 if(fenergy == 0. && elvlcstream.peek() !='\n')
00206 {
00207 elvlcstream.get(obtemp,lvl_nrg_col);
00208 if(atof(obtemp) != 0.)
00209 {
00210 fenergy = (realnum) atof(obtemp);
00211 }
00212
00213
00214 }
00215
00216 elvlcstream.ignore(INT_MAX,'\n');
00217 dBaseStatesEnergy[ipLev] = fenergy;
00218 }
00219 else
00220 {
00221 fprintf( ioQQQ, " The data file %s is corrupted .\n",chEnFilename);
00222 cdEXIT(EXIT_FAILURE);
00223 }
00224
00225
00226
00227 if(ipLev>0)
00228 {
00229 if (dBaseStatesEnergy[ipLev] < dBaseStatesEnergy[ipLev-1])
00230 {
00231 lgEneLevOrder = false;
00232 }
00233 }
00234 }
00235
00236
00237 intNewIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
00238
00239 if(lgEneLevOrder == false)
00240 {
00241
00242
00243 intOldIndex = (long int*)MALLOC((unsigned long)(nMolLevs)* sizeof(long int));
00244
00245 spsort(dBaseStatesEnergy,nMolLevs,intOldIndex,2,&interror);
00246
00247 for( long i=0; i<nMolLevs; i++ )
00248 {
00249 ASSERT( intOldIndex[i] < nMolLevs );
00250 intNewIndex[intOldIndex[i]] = i;
00251 }
00252 if(DEBUGSTATE)
00253 {
00254 for( long i=0; i<nMolLevs; i++ )
00255 {
00256 printf("The %ld value of the relation matrix is %ld \n",i,intNewIndex[i]);
00257 }
00258 for( long i=0; i<nMolLevs; i++ )
00259 {
00260 printf("The %ld value of the sorted energy array is %f \n",i,dBaseStatesEnergy[i]);
00261 }
00262 }
00263 free( intOldIndex );
00264 }
00265 else
00266 {
00267
00268 for( long i=0; i<nMolLevs; i++ )
00269 {
00270 intNewIndex[i] = i;
00271 }
00272 }
00273
00274
00275 for( long i=0; i<nMolLevs; i++ )
00276 {
00277 for( long j=i+1; j<nMolLevs; j++ )
00278 {
00279 ASSERT( intNewIndex[i] != intNewIndex[j] );
00280 }
00281 }
00282
00283
00284
00285 elvlcstream.seekg(0,ios::beg);
00286
00287
00288 const int lvl_skipto_statwt = 37;
00289
00290 const int lvl_statwt_col = 4;
00291
00292 for( long ipLevOld=0; ipLevOld<nMolLevs; ipLevOld++ )
00293 {
00294 long ipLevNew = intNewIndex[ipLevOld];
00295 char gtemp[lvl_statwt_col],thtemp[lvl_nrg_col],obtemp[lvl_nrg_col];
00296
00297
00298 strcpy(dBaseStates[intNS][ipLevNew].chLabel, " ");
00299 strncpy(dBaseStates[intNS][ipLevNew].chLabel,Species[intNS].chLabel, 4);
00300 dBaseStates[intNS][ipLevNew].chLabel[4] = '\0';
00301
00302
00303 dBaseStates[intNS][ipLevNew].sp = &Species[intNS];
00304
00305
00306 elvlcstream.seekg(lvl_skipto_statwt,ios::cur);
00307 elvlcstream.get(gtemp,lvl_statwt_col);
00308 fstatwt = (realnum)atof(gtemp);
00309
00310 if(fstatwt <= 0.)
00311 {
00312 fprintf( ioQQQ, " WARNING: A positive non zero value is expected for the statistical weight but was not found in %s\n", chEnFilename);
00313 fstatwt = 1.f;
00314
00315 }
00316 dBaseStates[intNS][ipLevNew].g = fstatwt;
00317
00318
00319 elvlcstream.get(thtemp,lvl_nrg_col);
00320 elvlcstream.seekg(lvl_skip_ryd,ios::cur);
00321 fenergy = (realnum) atof(thtemp);
00322
00323
00324
00325 if(fenergy == 0. && elvlcstream.peek() !='\n')
00326 {
00327 elvlcstream.get(obtemp,lvl_nrg_col);
00328 if(atof(obtemp) != 0.)
00329 {
00330 fenergy = (realnum) atof(obtemp);
00331 }
00332 }
00333 elvlcstream.ignore(INT_MAX,'\n');
00334 dBaseStates[intNS][ipLevNew].energy.set(fenergy,"cm^-1");
00335 }
00336
00337 elvlcstream.close();
00338
00339
00340 for( long ipHi=1; ipHi<nMolLevs; ipHi++ )
00341 {
00342 for( long ipLo=0; ipLo<ipHi; ipLo++ )
00343 {
00344 fenergyWN = (realnum)MAX2( 1E-10, dBaseStates[intNS][ipHi].energy.WN() - dBaseStates[intNS][ipLo].energy.WN() );
00345 fenergyK = (realnum)(fenergyWN*T1CM);
00346
00347 dBaseTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
00348 dBaseTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
00349 dBaseTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
00350 dBaseTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
00351 }
00352 }
00353
00354
00355
00356
00357
00358
00359 long wgfarows = -1;
00360 if (wgfastream.is_open())
00361 {
00362 int nj = 0;
00363 char otemp[eof_col];
00364 while(nj != -1)
00365 {
00366 wgfastream.get(otemp,eof_col);
00367 wgfastream.ignore(INT_MAX,'\n');
00368 nj = atoi(otemp);
00369 wgfarows++;
00370 }
00371 wgfastream.seekg(0,ios::beg);
00372 }
00373 else
00374 fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
00375
00376
00377 const int line_index_col = 6;
00378
00379 const int line_nrg_to_eina =15;
00380
00381 const int line_eina_col = 16;
00382 char lvltemp[line_index_col];
00383
00384 if (wgfastream.is_open())
00385 {
00386 for (long i = 0;i<wgfarows;i++)
00387 {
00388 wgfastream.get(lvltemp,line_index_col);
00389 long ipLo = atoi(lvltemp)-1;
00390 wgfastream.get(lvltemp,line_index_col);
00391 long ipHi = atoi(lvltemp)-1;
00392
00393 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
00394 {
00395
00396 wgfastream.ignore(INT_MAX,'\n');
00397 continue;
00398 }
00399
00400 if( ipHi == ipLo )
00401 {
00402 fprintf( ioQQQ," WARNING: Upper level = lower for a radiative transition in %s. Ignoring.\n", chTraFilename );
00403 wgfastream.ignore(INT_MAX,'\n');
00404 continue;
00405 }
00406
00407
00408 ipLo = intNewIndex[ipLo];
00409 ipHi = intNewIndex[ipHi];
00410
00411 ASSERT( ipHi != ipLo );
00412 ASSERT(ipHi >= 0);
00413 ASSERT(ipLo >= 0);
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426 if( ipHi < ipLo )
00427 {
00428 long ipLoTemp = ipLo;
00429 long ipHiTemp = ipHi;
00430 ipHi = max( ipHiTemp, ipLoTemp );
00431 ipLo = min( ipHiTemp, ipLoTemp );
00432
00433 }
00434
00435
00436 char trantemp[lvl_nrg_col];
00437 wgfastream.get(trantemp,lvl_nrg_col);
00438 fWLAng = (realnum)atof(trantemp);
00439
00440
00441
00442
00443 if( fWLAng <= 0. )
00444 {
00445
00446
00447 fWLAng = (realnum)(1e8/
00448 (dBaseStates[intNS][ipHi].energy.WN() - dBaseStates[intNS][ipLo].energy.WN()));
00449 }
00450
00451 wgfastream.seekg(line_nrg_to_eina,ios::cur);
00452 wgfastream.get(trantemp,line_eina_col);
00453 feinsteina = (realnum)atof(trantemp);
00454 if( feinsteina < 1e-20 )
00455 {
00456 static bool notPrintedYet = true;
00457 if( notPrintedYet )
00458 {
00459 fprintf( ioQQQ," WARNING: Radiative rate(s) equal to zero in %s.\n", chTraFilename );
00460 notPrintedYet = false;
00461 }
00462 wgfastream.ignore(INT_MAX,'\n');
00463 continue;
00464 }
00465
00466 fixit();
00467
00468 wgfastream.getline(chLine,INT_MAX);
00469 if( nMatch("auto", chLine) )
00470 {
00471 if( dBaseTrans[intNS][ipHi][ipLo].Emis != NULL )
00472 {
00473 dBaseTrans[intNS][ipHi][ipLo].Emis->AutoIonizFrac =
00474 feinsteina/(dBaseTrans[intNS][ipHi][ipLo].Emis->Aul + feinsteina);
00475 ASSERT( dBaseTrans[intNS][ipHi][ipLo].Emis->AutoIonizFrac >= 0. );
00476 ASSERT( dBaseTrans[intNS][ipHi][ipLo].Emis->AutoIonizFrac <= 1. );
00477 }
00478 continue;
00479 }
00480
00481 dBaseTrans[intNS][ipHi][ipLo].Emis
00482 = AddLine2Stack(feinsteina, &dBaseTrans[intNS][ipHi][ipLo]);
00483 dBaseTrans[intNS][ipHi][ipLo].WLAng = fWLAng;
00484
00485 fenergyWN = (realnum)(1e+8/fWLAng);
00486 fenergyK = (realnum)(fenergyWN*T1CM);
00487 dBaseTrans[intNS][ipHi][ipLo].EnergyK = fenergyK;
00488 dBaseTrans[intNS][ipHi][ipLo].EnergyWN = fenergyWN;
00489 dBaseTrans[intNS][ipHi][ipLo].EnergyErg = (realnum)ERG1CM *fenergyWN;
00490 dBaseTrans[intNS][ipHi][ipLo].WLAng = (realnum)(1e+8/fenergyWN/RefIndex(fenergyWN));
00491 dBaseTrans[intNS][ipHi][ipLo].Emis->gf = (realnum)GetGF(dBaseTrans[intNS][ipHi][ipLo].Emis->Aul,
00492 dBaseTrans[intNS][ipHi][ipLo].EnergyWN, dBaseTrans[intNS][ipHi][ipLo].Hi->g);
00493
00494 if(DEBUGSTATE)
00495 {
00496 fprintf( ioQQQ, "The lower level is %ld \n",ipLo);
00497 fprintf( ioQQQ, "The upper level is %ld \n",ipHi);
00498 fprintf( ioQQQ, "The Einstein A is %f \n",dBaseTrans[intNS][ipHi][ipLo].Emis->Aul);
00499 fprintf( ioQQQ, "The wavelength of the transition is %f \n",dBaseTrans[intNS][ipHi][ipLo].WLAng );
00500 }
00501 }
00502 }
00503 else fprintf( ioQQQ, " The data file %s is corrupted .\n",chTraFilename);
00504 wgfastream.close();
00505
00506
00507 AtmolCollSplines[intNS] = (CollSplinesArray***)MALLOC(nMolLevs *sizeof(CollSplinesArray**));
00508 for( long ipHi=0; ipHi<nMolLevs; ipHi++ )
00509 {
00510 AtmolCollSplines[intNS][ipHi] = (CollSplinesArray **)MALLOC((unsigned long)(nMolLevs)*sizeof(CollSplinesArray *));
00511 for( long ipLo=0; ipLo<nMolLevs; ipLo++ )
00512 {
00513 AtmolCollSplines[intNS][ipHi][ipLo] =
00514 (CollSplinesArray *)MALLOC((unsigned long)(NUM_COLLIDERS)*sizeof(CollSplinesArray ));
00515
00516 for( long k=0; k<NUM_COLLIDERS; k++ )
00517 {
00518
00519 AtmolCollSplines[intNS][ipHi][ipLo][k].collspline = NULL;
00520 AtmolCollSplines[intNS][ipHi][ipLo][k].SplineSecDer = NULL;
00521 AtmolCollSplines[intNS][ipHi][ipLo][k].nSplinePts = -1;
00522 AtmolCollSplines[intNS][ipHi][ipLo][k].intTranType = -1;
00523 AtmolCollSplines[intNS][ipHi][ipLo][k].gf = BIGDOUBLE;
00524 AtmolCollSplines[intNS][ipHi][ipLo][k].EnergyDiff = BIGDOUBLE;
00525 AtmolCollSplines[intNS][ipHi][ipLo][k].ScalingParam = BIGDOUBLE;
00526 }
00527 }
00528 }
00529
00530
00531
00532
00533
00534
00535 for( long ipCollider=0; ipCollider<=1; ipCollider++ )
00536 {
00537 char chFilename[FILENAME_PATH_LENGTH_2];
00538
00539 if( ipCollider==0 )
00540 {
00541 strcpy( chFilename, chEleColFilename );
00542 }
00543 else if( ipCollider==1 )
00544 {
00545 if( !lgProtonData )
00546 break;
00547 fprintf( ioQQQ," WARNING: Chianti proton collision data have different format than electron data files!\n" );
00548 strcpy( chFilename, chProColFilename );
00549 }
00550 else
00551 TotalInsanity();
00552
00553
00554 strcpy(chLine,"A");
00555
00556 fstream splupsstream;
00557 open_data(splupsstream, chFilename,mode_r);
00558
00559
00560 const int cs_eof_col = 4;
00561
00562 const int cs_index_col = 4;
00563
00564 const int cs_trantype_col = 4;
00565
00566
00567 const int cs_values_col = 11;
00568
00569 if (splupsstream.is_open())
00570 {
00571 int nj = 0;
00572
00573 long splupslines = -1;
00574 char otemp[cs_eof_col];
00575 while(nj != -1)
00576 {
00577 splupsstream.get(otemp,cs_eof_col);
00578 splupsstream.ignore(INT_MAX,'\n');
00579 nj = atoi(otemp);
00580 splupslines++;
00581 }
00582 splupsstream.seekg(0,ios::beg);
00583
00584 for (int m = 0;m<splupslines;m++)
00585 {
00586 if( ipCollider == 0 )
00587 {
00588 splupsstream.seekg(6,ios::cur);
00589 }
00590
00591
00592 splupsstream.get(otemp,cs_index_col);
00593 long ipLo = atoi(otemp)-1;
00594 splupsstream.get(otemp,cs_index_col);
00595 long ipHi = atoi(otemp)-1;
00596
00597 if( ipLo >= nMolLevs || ipHi >= nMolLevs )
00598 {
00599
00600 splupsstream.ignore(INT_MAX,'\n');
00601 continue;
00602 }
00603
00604 if( ipLo >= HighestIndexInFile || ipHi >= HighestIndexInFile )
00605 {
00606 fprintf( ioQQQ," WARNING: Problem with indices in CHIANTI file %s. Often due to indices > 100 in fixed format. Disabling model.\n", chFilename );
00607 fixit();
00608 break;
00609 }
00610
00611 ASSERT( ipLo < nMolLevs );
00612 ASSERT( ipHi < nMolLevs );
00613
00614 splupsstream.get(otemp,cs_trantype_col);
00615 intTranType = atoi(otemp);
00616
00617 char qtemp[cs_values_col];
00618 splupsstream.get(qtemp,cs_values_col);
00619 fGF = atof(qtemp);
00620
00621 splupsstream.get(qtemp,cs_values_col);
00622 fEnergyDiff = atof(qtemp);
00623
00624 splupsstream.get(qtemp,cs_values_col);
00625 fScalingParam = atof(qtemp);
00626
00627
00628 ipLo = intNewIndex[ipLo];
00629 ipHi = intNewIndex[ipHi];
00630
00631 ASSERT( ipLo < nMolLevs );
00632 ASSERT( ipHi < nMolLevs );
00633 ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline == NULL );
00634 ASSERT( AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer == NULL );
00635
00636 const int CHIANTI_SPLINE_MAX=9, CHIANTI_SPLINE_MIN=5;
00637 STATIC_ASSERT(CHIANTI_SPLINE_MAX > CHIANTI_SPLINE_MIN);
00638
00639
00640 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline =
00641 (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
00642 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer =
00643 (double *)MALLOC((unsigned long)(CHIANTI_SPLINE_MAX)*sizeof(double));
00644
00645
00646 for( intsplinepts=0; intsplinepts<=CHIANTI_SPLINE_MAX; intsplinepts++ )
00647 {
00648
00649 char p = splupsstream.peek();
00650 if( p == '\n' )
00651 {
00652
00653
00654 if( (intsplinepts != CHIANTI_SPLINE_MAX) && (intsplinepts != CHIANTI_SPLINE_MIN) )
00655 {
00656 static bool notPrintedYet = true;
00657 if( notPrintedYet )
00658 {
00659 fprintf( ioQQQ, " WARNING: Wrong number of spline points in %s.\n", chFilename);
00660 notPrintedYet = false;
00661 }
00662 for( long j=intsplinepts-1; j<CHIANTI_SPLINE_MAX; j++ )
00663 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[j] = 0.;
00664 }
00665 else
00666 break;
00667 }
00668 else
00669 {
00670 if( intsplinepts >= CHIANTI_SPLINE_MAX )
00671 {
00672 fprintf( ioQQQ, " WARNING: More spline points than expected in %s, indices %3li %3li. Ignoring extras.\n", chFilename, ipHi, ipLo );
00673 break;
00674 }
00675 ASSERT( intsplinepts < CHIANTI_SPLINE_MAX );
00676 double temp;
00677
00678 splupsstream.get(qtemp,cs_values_col);
00679 temp = atof(qtemp);
00680 if(temp < 0)
00681 {
00682 temp = 0.;
00683 }
00684 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intsplinepts] = temp;
00685 }
00686
00687 }
00688
00689 ASSERT( (intsplinepts == CHIANTI_SPLINE_MAX) || (intsplinepts == CHIANTI_SPLINE_MIN) );
00690
00691
00692 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].nSplinePts = intsplinepts;
00693
00694 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].intTranType = intTranType;
00695
00696 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].gf = fGF;
00697
00698 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].EnergyDiff = fEnergyDiff;
00699
00700 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].ScalingParam = fScalingParam;
00701
00702
00703
00704 xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
00705 spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
00706 spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
00707
00708
00709 if(intsplinepts == CHIANTI_SPLINE_MIN)
00710 {
00711 for(intxs=0;intxs<CHIANTI_SPLINE_MIN;intxs++)
00712 {
00713 xs[intxs] = 0.25*intxs;
00714 spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
00715 }
00716 }
00717 else if(intsplinepts == CHIANTI_SPLINE_MAX)
00718 {
00719 for(intxs=0;intxs<CHIANTI_SPLINE_MAX;intxs++)
00720 {
00721 xs[intxs] = 0.125*intxs;
00722 spl[intxs] = AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].collspline[intxs];
00723 }
00724 }
00725 else
00726 TotalInsanity();
00727
00728 spline(xs, spl,intsplinepts,2e31,2e31,spl2);
00729
00730
00731 for( long i=0; i<intsplinepts; i++)
00732 {
00733 AtmolCollSplines[intNS][ipHi][ipLo][ipCollider].SplineSecDer[i] = spl2[i];
00734 }
00735
00736 free( xs );
00737 free( spl );
00738 free( spl2 );
00739 splupsstream.ignore(INT_MAX,'\n');
00740 }
00741 splupsstream.close();
00742 }
00743
00744 }
00745
00746
00747 free( dBaseStatesEnergy );
00748 free( intNewIndex );
00749
00750
00751 fclose( ioElecCollData );
00752 if( lgProtonData )
00753 fclose( ioProtCollData );
00754
00755 return;
00756 }