00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "cddrive.h"
00007 #include "physconst.h"
00008 #include "geometry.h"
00009 #include "radius.h"
00010 #include "rfield.h"
00011 #include "opacity.h"
00012 #include "optimize.h"
00013 #include "grid.h"
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 static long int iplo , iphi;
00025
00026
00027 static double Elo , Ehi;
00028
00029
00030 static double Spec_cont( float cont[] );
00031
00032 void cdSPEC(
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050 int nOption ,
00051
00052
00053
00054 double EnergyLow[] ,
00055
00056
00057 long int nEnergy ,
00058
00059
00060 double ReturnedSpectrum[] )
00061
00062 {
00063
00064 float *cont ,
00065 refac;
00066 long int ncell , j;
00067
00068
00069 bool lgFREE;
00070
00071 DEBUG_ENTRY( "cdSPEC()" );
00072
00073 if( nOption == 1 )
00074 {
00075
00076 cont = rfield.FluxSave;
00077 lgFREE = false;
00078 }
00079 else if( nOption == 2 )
00080 {
00081
00082
00083 cont = rfield.flux;
00084 lgFREE = false;
00085 }
00086 else if( nOption == 3 )
00087 {
00088
00089 lgFREE = false;
00090 cont = rfield.ConRefIncid;
00091 }
00092 else if( nOption == 4 )
00093 {
00094
00095 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00096
00097
00098 lgFREE = true;
00099 refac = (float)radius.r1r0sq*geometry.covgeo;
00100 for( j=0; j<rfield.nflux; ++j)
00101 {
00102 cont[j] = rfield.ConEmitOut[j]*refac;
00103 }
00104 }
00105 else if( nOption == 5 )
00106 {
00107
00108 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00109
00110 lgFREE = true;
00111 refac = (float)radius.r1r0sq*geometry.covgeo;
00112 for( j=0; j<rfield.nflux; ++j)
00113 {
00114 cont[j] = rfield.ConEmitReflec[j]*refac;
00115 }
00116 }
00117 else if( nOption == 6 )
00118 {
00119
00120 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00121
00122 lgFREE = true;
00123
00124 refac = (float)radius.r1r0sq*geometry.covgeo;
00125 for( j=0; j<rfield.nflux; ++j)
00126 {
00127
00128
00129 cont[j] = rfield.outlin[j] *rfield.widflx[j]/rfield.anu[j]*refac;
00130 }
00131 }
00132 else if( nOption == 7 )
00133 {
00134
00135 if( geometry.lgSphere )
00136 {
00137 refac = 0.;
00138 }
00139 else
00140 {
00141 refac = 1.;
00142 }
00143
00144 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00145
00146 lgFREE = true;
00147 for( j=0; j<rfield.nflux; ++j)
00148 {
00149
00150
00151 cont[j] = rfield.reflin[j] *rfield.widflx[j]/rfield.anu[j]*refac;
00152 }
00153 }
00154 else
00155 {
00156 fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
00157 puts( "[Stop in cdSPEC]" );
00158 cdEXIT(EXIT_FAILURE);
00159 }
00160
00161
00162 iplo = 0;
00163 iphi = 0;
00164
00165 for( ncell = 0; ncell < nEnergy-1; ++ncell )
00166 {
00167
00168 Elo = EnergyLow[ncell];
00169 Ehi = EnergyLow[ncell+1];
00170 ReturnedSpectrum[ncell] = Spec_cont( cont );
00171 }
00172
00173
00174 if( lgFREE )
00175 {
00176 free(cont);
00177 }
00178 DEBUG_EXIT( "cdSPEC()" );
00179 return;
00180 }
00181
00182
00183 static double Spec_cont( float cont[] )
00184 {
00185 double sum;
00186 long int i;
00187
00188
00189
00190
00191
00192 if( iplo >= rfield.nflux )
00193 return 0.;
00194
00195
00196
00197
00198 while( (iplo < rfield.nflux) &&
00199 !(Elo <= (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) &&
00200 (Elo >= (rfield.AnuOrg[iplo]-rfield.widflx[iplo]/2.) )) )
00201 ++iplo;
00202
00203 iphi = iplo;
00204
00205 while( (iphi < rfield.nflux) &&
00206 !(Ehi <= (rfield.AnuOrg[iphi+1]-rfield.widflx[iphi+1]/2.) &&
00207 (Ehi >= (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) )) )
00208 ++iphi;
00209
00210 sum = 0.;
00211 if( iphi-iplo>1 )
00212 {
00213
00214
00215 for( i=iplo+1; i<iphi; ++i )
00216 {
00217 sum += cont[i] * rfield.anu2[i];
00218 }
00219 }
00220 ASSERT( sum>=0.);
00221
00222
00223
00224
00225 if( iplo == iphi )
00226 {
00227
00228 sum = cont[iplo]/rfield.widflx[iplo]*(Ehi-Elo) * rfield.anu2[iplo];
00229 ASSERT( sum>=0.);
00230 }
00231 else
00232 {
00233
00234 double piece;
00235
00236
00237 piece = cont[iplo]/rfield.widflx[iplo]*( (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) - Elo ) *
00238 rfield.anu2[iplo];
00239 ASSERT( piece >= 0.);
00240 sum += piece;
00241
00242
00243 piece = cont[iphi]/rfield.widflx[iphi]*(Ehi - (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) ) *
00244 rfield.anu2[iphi];
00245 ASSERT( piece >= 0.);
00246 sum += piece;
00247 ASSERT( sum>=0.);
00248 }
00249
00250
00251 sum *= EN1RYD / (Ehi - Elo );
00252
00253 ASSERT( sum >=0. );
00254 return sum;
00255
00256 # if 0
00257
00258
00259
00260
00261 if( (Ehi-Elo) < rfield.widflx[ip] )
00262 {
00263
00264 double Emiddle = (Elo + Ehi)/2.;
00265 double f1 , f2 ,finterp;
00266 ASSERT( (Elo >= (rfield.anu[ip]+rfield.widflx[ip]/2.
00267
00268
00269
00270
00271 f1 = cont[ip-1] /rfield.widflx[ip-1] *rfield.anu2[ip-1]*EN1RYD;
00272 f2 = cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD;
00273
00274 finterp = f1 + (f2-f1) /rfield.widflx[ip] *
00275 fabs(Emiddle-(rfield.anu[ip]-rfield.widflx[ip]/2.));
00276 return( finterp );
00277 }
00278 else
00279 {
00280
00281 sum = cont[ip] * rfield.anu2[ip]* EN1RYD;
00282
00283 EcHi = ( rfield.anu[ip] + rfield.widflx[ip]/2. + rfield.anu[ip+1] - rfield.widflx[ip+1]/2.)/2.;
00284
00285 sum *= (EcHi - Elo )/rfield.widflx[ip];
00286 }
00287
00288
00289 ++ip;
00290
00291
00292
00293 while( Ehi > rfield.anu[ip]+rfield.widflx[ip]/2. )
00294 {
00295 sum += cont[ip] * rfield.anu2[ip]*EN1RYD;
00296 ++ip;
00297 }
00298
00299
00300 EcLo = ( rfield.anu[ip] - rfield.widflx[ip]/2. + rfield.anu[ip-1] + rfield.widflx[ip-1]/2.)/2.;
00301
00302
00303 sum += cont[ip] * rfield.anu2[ip] * EN1RYD *
00304 (Ehi - EcLo )/rfield.widflx[ip];
00305
00306
00307 sum /= (Ehi - Elo );
00308
00309 return(sum);
00310 # endif
00311 }
00312
00313
00314 void cdSPEC2(
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 int nOption ,
00335
00336
00337 long int nEnergy ,
00338
00339
00340 float ReturnedSpectrum[] )
00341
00342 {
00343
00344 float *cont ,
00345 refac;
00346 long int ncell , j;
00347
00348
00349 bool lgFREE;
00350
00351 DEBUG_ENTRY( "cdSPEC2()" );
00352
00353 ASSERT( nOption <= NUM_OUTPUT_TYPES );
00354
00355 if( nOption == 1 )
00356 {
00357
00358 cont = rfield.FluxSave;
00359 lgFREE = false;
00360 }
00361 else if( nOption == 2 )
00362 {
00363
00364
00365 cont = rfield.flux;
00366 lgFREE = false;
00367 }
00368 else if( nOption == 3 )
00369 {
00370
00371 lgFREE = false;
00372 cont = rfield.ConRefIncid;
00373 }
00374 else if( nOption == 4 )
00375 {
00376
00377 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00378
00379 lgFREE = true;
00380 refac = (float)radius.r1r0sq*geometry.covgeo;
00381 for( j=0; j<rfield.nflux; ++j)
00382 {
00383 cont[j] = rfield.ConEmitOut[j]*refac;
00384 }
00385 }
00386 else if( nOption == 5 )
00387 {
00388
00389 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00390
00391 lgFREE = true;
00392 for( j=0; j<rfield.nflux; ++j)
00393 {
00394 cont[j] = rfield.ConEmitReflec[j];
00395 }
00396 }
00397 else if( nOption == 6 )
00398 {
00399
00400 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00401
00402 lgFREE = true;
00403
00404 refac = (float)radius.r1r0sq*geometry.covgeo;
00405 for( j=0; j<rfield.nflux; ++j)
00406 {
00407 cont[j] = rfield.outlin[j]*refac;
00408 }
00409 }
00410 else if( nOption == 7 )
00411 {
00412
00413 if( geometry.lgSphere )
00414 {
00415 refac = 0.;
00416 }
00417 else
00418 {
00419 refac = 1.;
00420 }
00421
00422 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00423
00424 lgFREE = true;
00425 for( j=0; j<rfield.nflux; ++j)
00426 {
00427
00428
00429 cont[j] = rfield.reflin[j]*refac;
00430 }
00431 }
00432 else if( nOption == 8 )
00433 {
00434
00435 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00436
00437 lgFREE = true;
00438
00439 refac = (float)radius.r1r0sq*geometry.covgeo;
00440 for( j=0; j<rfield.nflux; ++j)
00441 {
00442
00443
00444 cont[j] = (rfield.ConEmitOut[j]+ rfield.outlin[j])*refac
00445 + rfield.flux[j]*(float)radius.r1r0sq;
00446 }
00447 }
00448 else if( nOption == 9 )
00449 {
00450
00451 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00452
00453 lgFREE = true;
00454 for( j=0; j<rfield.nflux; ++j)
00455 {
00456
00457
00458 cont[j] = ((rfield.ConRefIncid[j]+rfield.ConEmitReflec[j]) +
00459 rfield.reflin[j]);
00460 }
00461 }
00462 else if( nOption == 10 )
00463 {
00464
00465 cont = (float*)MALLOC( sizeof(float)*(size_t)rfield.nupper );
00466
00467 lgFREE = true;
00468 for( j=0; j<nEnergy; ++j)
00469 {
00470
00471
00472 cont[j] = opac.ExpmTau[j]*rfield.trans_coef_total[j];
00473 }
00474 }
00475 else
00476 {
00477 fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
00478 puts( "[Stop in cdSPEC]" );
00479 cdEXIT(EXIT_FAILURE);
00480 }
00481
00482
00483 iplo = 0;
00484 iphi = 0;
00485
00486 for( ncell = 0; ncell < nEnergy; ++ncell )
00487 {
00488 if( ncell >= rfield.nflux )
00489 {
00490 ReturnedSpectrum[ncell] = SMALLFLOAT;
00491 }
00492 else
00493 {
00494 ReturnedSpectrum[ncell] = cont[ncell];
00495 }
00496 ASSERT( ReturnedSpectrum[ncell] >=0.f );
00497 }
00498
00499
00500 if( lgFREE )
00501 {
00502 free(cont);
00503 }
00504 DEBUG_EXIT( "cdSPEC2()" );
00505 return;
00506 }