cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*cdSPEC returns the spectrum needed for Keith Arnaud's XSPEC */ 00004 /*Spec_cont called by cdSPEC to generate actual spectrum */ 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 "grid.h" 00013 00014 /* 00015 * this routine returns the spectrum needed for Keith Arnaud's XSPEC 00016 * X-Ray analysis code. It should be called after cdDrive has successfully 00017 * computed a model. the calling routine must ensure that the vectors 00018 * have enough space to store the resulting spectrum, 00019 * given the bounds and energy resolution 00020 */ 00021 00022 /* array index within energy array in Cloudy grids - this has file scope */ 00023 static long int iplo , iphi; 00024 00025 /* energies of lower and upper bound of XSPEC cell we are considering */ 00026 static double Elo , Ehi; 00027 00028 /* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */ 00029 STATIC double Spec_cont( realnum cont[] ); 00030 00031 void cdSPEC( 00032 /* option - the type of spectrum to be returned 00033 * 1 the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1 00034 * 00035 * 2 the attenuated incident continuum, same units 00036 * 3 the reflected continuum, same units 00037 * 00038 * 4 diffuse continuous emission outward direction 00039 * 5 diffuse continuous emission, reflected 00040 * 00041 * 6 collisional+recombination lines, outward 00042 * 7 collisional+recombination lines, reflected 00043 * 00044 * 8 pumped lines, outward <= not implemented 00045 * 9 pumped lines, reflected <= not implemented 00046 * 00047 * all lines and continuum emitted by the cloud assume full coverage of 00048 * continuum source */ 00049 int nOption , 00050 00051 /* the energy of the lower edge of each cell 00052 * (in Ryd to be consistent with all the rest of Cloudy) */ 00053 double EnergyLow[] , 00054 00055 /* the number of cells + 1*/ 00056 long int nEnergy , 00057 00058 /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */ 00059 double ReturnedSpectrum[] ) 00060 00061 { 00062 /* this pointer will bet set to one of the cloudy continuum arrays */ 00063 realnum *cont , 00064 refac; 00065 long int ncell , j; 00066 00067 /* flag saying whether we will need to free cont at the end */ 00068 bool lgFREE; 00069 00070 DEBUG_ENTRY( "cdSPEC()" ); 00071 00072 if( nOption == 1 ) 00073 { 00074 /* this is the incident continuum, col 2 of punch continuum command */ 00075 cont = rfield.flux_total_incident; 00076 lgFREE = false; 00077 } 00078 else if( nOption == 2 ) 00079 { 00080 /* the attenuated transmitted continuum, no diffuse emission, 00081 * col 3 of punch continuum command */ 00082 cont = rfield.flux; 00083 lgFREE = false; 00084 } 00085 else if( nOption == 3 ) 00086 { 00087 /* reflected incident continuum, col 6 of punch continuum command */ 00088 lgFREE = false; 00089 cont = rfield.ConRefIncid; 00090 } 00091 else if( nOption == 4 ) 00092 { 00093 /* diffuse continuous emission outward direction */ 00094 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00095 00096 /* need to free the vector once done */ 00097 lgFREE = true; 00098 refac = (realnum)radius.r1r0sq*geometry.covgeo; 00099 for( j=0; j<rfield.nflux; ++j) 00100 { 00101 cont[j] = rfield.ConEmitOut[j]*refac; 00102 } 00103 } 00104 else if( nOption == 5 ) 00105 { 00106 /* reflected diffuse continuous emission */ 00107 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00108 /* need to free the vector once done */ 00109 lgFREE = true; 00110 refac = (realnum)radius.r1r0sq*geometry.covgeo; 00111 for( j=0; j<rfield.nflux; ++j) 00112 { 00113 cont[j] = rfield.ConEmitReflec[j]*refac; 00114 } 00115 } 00116 else if( nOption == 6 ) 00117 { 00118 /* all outward lines, */ 00119 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00120 /* need to free the vector once done */ 00121 lgFREE = true; 00122 /* correct back to inner radius */ 00123 refac = (realnum)radius.r1r0sq*geometry.covgeo; 00124 for( j=0; j<rfield.nflux; ++j) 00125 { 00126 /* units of lines here are to cancel out with tricks applied to continuum cells 00127 * when finally extracted below */ 00128 cont[j] = rfield.outlin[j] *rfield.widflx[j]/rfield.anu[j]*refac; 00129 } 00130 } 00131 else if( nOption == 7 ) 00132 { 00133 /* all reflected lines */ 00134 if( geometry.lgSphere ) 00135 { 00136 refac = 0.; 00137 } 00138 else 00139 { 00140 refac = 1.; 00141 } 00142 00143 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00144 /* need to free the vector once done */ 00145 lgFREE = true; 00146 for( j=0; j<rfield.nflux; ++j) 00147 { 00148 /* units of lines here are to cancel out with tricks applied to continuum cells 00149 * when finally extracted below */ 00150 cont[j] = rfield.reflin[j] *rfield.widflx[j]/rfield.anu[j]*refac; 00151 } 00152 } 00153 else 00154 { 00155 fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption); 00156 cdEXIT(EXIT_FAILURE); 00157 } 00158 00159 /* this is array index used in Spec_cont */ 00160 iplo = 0; 00161 iphi = 0; 00162 /* now generate the continua */ 00163 for( ncell = 0; ncell < nEnergy-1; ++ncell ) 00164 { 00165 /* Spec_cont will find the spectrum between these two energies */ 00166 Elo = EnergyLow[ncell]; 00167 Ehi = EnergyLow[ncell+1]; 00168 ReturnedSpectrum[ncell] = Spec_cont( cont ); 00169 } 00170 00171 /* need to free the vector once done if this flag was set */ 00172 if( lgFREE ) 00173 { 00174 free(cont); 00175 } 00176 return; 00177 } 00178 00179 /* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */ 00180 STATIC double Spec_cont( realnum cont[] ) 00181 { 00182 double sum; 00183 long int i; 00184 00185 DEBUG_ENTRY( "Spec_cont()" ); 00186 00187 /* iplo and iphi are set to zero before this is called to generate a full spectrum, 00188 * since energies are in increasing order, after the first call it should only 00189 * be necessary to increase ip above where it now is */ 00190 /* return zero if continuum does not go this high */ 00191 if( iplo >= rfield.nflux ) 00192 return 0.; 00193 00194 /* now skip out to where continuum starts, often we will already be there 00195 * if this is a successive call */ 00196 /* iplo should be index where Elo is greater than anu-widflx/1. */ 00197 while( (iplo < rfield.nflux) && 00198 !(Elo <= (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) && 00199 (Elo >= (rfield.AnuOrg[iplo]-rfield.widflx[iplo]/2.) )) ) 00200 ++iplo; 00201 00202 iphi = iplo; 00203 /* iphi should be index where Ehi is greater than anu-widflx/1. */ 00204 while( (iphi < rfield.nflux) && 00205 !(Ehi <= (rfield.AnuOrg[iphi+1]-rfield.widflx[iphi+1]/2.) && 00206 (Ehi >= (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) )) ) 00207 ++iphi; 00208 00209 sum = 0.; 00210 if( iphi-iplo>1 ) 00211 { 00212 /* this branch when more than two cells are involved - 00213 * begin by summing intermediate cells */ 00214 for( i=iplo+1; i<iphi; ++i ) 00215 { 00216 sum += cont[i] * rfield.anu2[i]; 00217 } 00218 } 00219 ASSERT( sum>=0.); 00220 00221 /* at this point sum will often be zero, and we need to add on flux from 00222 * the iplo and iphi cells - there are two possibilities - iplo and iphi can 00223 * be the same cell, or different cells */ 00224 if( iplo == iphi ) 00225 { 00226 /* we want only a piece of the ip cell */ 00227 sum = cont[iplo]/rfield.widflx[iplo]*(Ehi-Elo) * rfield.anu2[iplo]; 00228 ASSERT( sum>=0.); 00229 } 00230 else 00231 { 00232 /* we want to add on a fraction of the low and high energy cells */ 00233 double piece; 00234 00235 /* the low energy cell */ 00236 piece = cont[iplo]/rfield.widflx[iplo]*( (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) - Elo ) * 00237 rfield.anu2[iplo]; 00238 ASSERT( piece >= 0.); 00239 sum += piece; 00240 00241 /* now the high energy cell */ 00242 piece = cont[iphi]/rfield.widflx[iphi]*(Ehi - (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) ) * 00243 rfield.anu2[iphi]; 00244 ASSERT( piece >= 0.); 00245 sum += piece; 00246 ASSERT( sum>=0.); 00247 } 00248 00249 /* now divide by total energy width */ 00250 sum *= EN1RYD / (Ehi - Elo ); 00251 00252 ASSERT( sum >=0. ); 00253 return sum; 00254 00255 # if 0 00256 /* at this point Elo is less than the upper edge energy of the ip cell */ 00257 /* >>chng 01 jul 23, change logic, use this branch if desired width smaller than 00258 * intrinsic cloudy width */ 00259 /*if( rfield.anu[ip]+rfield.widflx[ip]/2. > Ehi )*/ 00260 if( (Ehi-Elo) < rfield.widflx[ip] ) 00261 { 00262 /* this is where we want the flux evaluated */ 00263 double Emiddle = (Elo + Ehi)/2.; 00264 double f1 , f2 ,finterp; 00265 ASSERT( (Elo >= (rfield.anu[ip]+rfield.widflx[ip]/2. 00266 /* this is case where requested cell is within a single Cloudy cell - 00267 * do linear interpolation */ 00268 /* >>chng 01 jul 23, change to linear interpolation */ 00269 /*return( cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD );*/ 00270 f1 = cont[ip-1] /rfield.widflx[ip-1] *rfield.anu2[ip-1]*EN1RYD; 00271 f2 = cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD; 00272 00273 finterp = f1 + (f2-f1) /rfield.widflx[ip] * 00274 fabs(Emiddle-(rfield.anu[ip]-rfield.widflx[ip]/2.)); 00275 return( finterp ); 00276 } 00277 else 00278 { 00279 /* this is case where we will add on more than one Cloudy cell */ 00280 sum = cont[ip] * rfield.anu2[ip]* EN1RYD; 00281 /* energy of upper edge of Cloudy cell */ 00282 EcHi = ( rfield.anu[ip] + rfield.widflx[ip]/2. + rfield.anu[ip+1] - rfield.widflx[ip+1]/2.)/2.; 00283 /* mutiply this by energy width of this first cell that counts to final integral */ 00284 sum *= (EcHi - Elo )/rfield.widflx[ip]; 00285 } 00286 00287 /* increment index since we will now work on next cell */ 00288 ++ip; 00289 00290 /* now add up total cells, this loop for cloudy cells totally within desired bins, 00291 * we will divide by energy width later, so factor of widflx is missing here */ 00292 while( Ehi > rfield.anu[ip]+rfield.widflx[ip]/2. ) 00293 { 00294 sum += cont[ip] * rfield.anu2[ip]*EN1RYD; 00295 ++ip; 00296 } 00297 00298 /* energy of upper edge of Cloudy cell */ 00299 EcLo = ( rfield.anu[ip] - rfield.widflx[ip]/2. + rfield.anu[ip-1] + rfield.widflx[ip-1]/2.)/2.; 00300 00301 /* now finish up with last cell, whose upper bound is higher than desired bin */ 00302 sum += cont[ip] * rfield.anu2[ip] * EN1RYD * 00303 (Ehi - EcLo )/rfield.widflx[ip]; 00304 00305 /* now divide by total energy width */ 00306 sum /= (Ehi - Elo ); 00307 00308 return(sum); 00309 # endif 00310 } 00311 00312 /* returns in units photons/cm^2/s/bin */ 00313 void cdSPEC2( 00314 /* option - the type of spectrum to be returned 00315 * 1 the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1 00316 * 00317 * 2 the attenuated incident continuum, same units 00318 * 3 the reflected continuum, same units 00319 * 00320 * 4 diffuse continuous emission outward direction 00321 * 5 diffuse continuous emission, reflected 00322 * 00323 * 6 collisional+recombination lines, outward 00324 * 7 collisional+recombination lines, reflected 00325 * 00326 * 8 total transmitted, lines and continuum 00327 * 9 total reflected, lines and continuum 00328 * 00329 *10 exp(-tau) to the illuminated face 00330 * 00331 * all lines and continuum emitted by the cloud assume full coverage of 00332 * continuum source */ 00333 int nOption , 00334 00335 /* the number of cells + 1*/ 00336 long int nEnergy , 00337 00338 /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */ 00339 realnum ReturnedSpectrum[] ) 00340 00341 { 00342 /* this pointer will bet set to one of the cloudy continuum arrays */ 00343 realnum *cont , 00344 refac; 00345 long int ncell , j; 00346 00347 /* flag saying whether we will need to free cont at the end */ 00348 bool lgFREE; 00349 00350 DEBUG_ENTRY( "cdSPEC2()" ); 00351 00352 ASSERT( nOption <= NUM_OUTPUT_TYPES ); 00353 00354 if( nOption == 1 ) 00355 { 00356 /* this is the incident continuum, col 2 of punch continuum command */ 00357 cont = rfield.flux_total_incident; 00358 lgFREE = false; 00359 } 00360 else if( nOption == 2 ) 00361 { 00362 /* the attenuated transmitted continuum, no diffuse emission, 00363 * col 3 of punch continuum command */ 00364 cont = rfield.flux; 00365 lgFREE = false; 00366 } 00367 else if( nOption == 3 ) 00368 { 00369 /* reflected incident continuum, col 6 of punch continuum command */ 00370 lgFREE = false; 00371 cont = rfield.ConRefIncid; 00372 } 00373 else if( nOption == 4 ) 00374 { 00375 /* diffuse continuous emission outward direction */ 00376 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00377 /* need to free the vector once done */ 00378 lgFREE = true; 00379 refac = (realnum)radius.r1r0sq*geometry.covgeo; 00380 for( j=0; j<rfield.nflux; ++j) 00381 { 00382 cont[j] = rfield.ConEmitOut[j]*refac; 00383 } 00384 } 00385 else if( nOption == 5 ) 00386 { 00387 /* reflected diffuse continuous emission */ 00388 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00389 /* need to free the vector once done */ 00390 lgFREE = true; 00391 for( j=0; j<rfield.nflux; ++j) 00392 { 00393 cont[j] = rfield.ConEmitReflec[j]; 00394 } 00395 } 00396 else if( nOption == 6 ) 00397 { 00398 /* all outward lines, */ 00399 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00400 /* need to free the vector once done */ 00401 lgFREE = true; 00402 /* correct back to inner radius */ 00403 refac = (realnum)radius.r1r0sq*geometry.covgeo; 00404 for( j=0; j<rfield.nflux; ++j) 00405 { 00406 cont[j] = rfield.outlin[j]*refac; 00407 } 00408 } 00409 else if( nOption == 7 ) 00410 { 00411 /* all reflected lines */ 00412 if( geometry.lgSphere ) 00413 { 00414 refac = 0.; 00415 } 00416 else 00417 { 00418 refac = 1.; 00419 } 00420 00421 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00422 /* need to free the vector once done */ 00423 lgFREE = true; 00424 for( j=0; j<rfield.nflux; ++j) 00425 { 00426 /* units of lines here are to cancel out with tricks applied to continuum cells 00427 * when finally extracted below */ 00428 cont[j] = rfield.reflin[j]*refac; 00429 } 00430 } 00431 else if( nOption == 8 ) 00432 { 00433 /* total transmitted continuum */ 00434 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00435 /* need to free the vector once done */ 00436 lgFREE = true; 00437 /* correct back to inner radius */ 00438 refac = (realnum)radius.r1r0sq*geometry.covgeo; 00439 for( j=0; j<rfield.nflux; ++j) 00440 { 00441 /* units of lines here are to cancel out with tricks applied to continuum cells 00442 * when finally extracted below */ 00443 cont[j] = (rfield.ConEmitOut[j]+ rfield.outlin[j])*refac 00444 + rfield.flux[j]*(realnum)radius.r1r0sq; 00445 } 00446 } 00447 else if( nOption == 9 ) 00448 { 00449 /* total reflected continuum */ 00450 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00451 /* need to free the vector once done */ 00452 lgFREE = true; 00453 for( j=0; j<rfield.nflux; ++j) 00454 { 00455 /* units of lines here are to cancel out with tricks applied to continuum cells 00456 * when finally extracted below */ 00457 cont[j] = ((rfield.ConRefIncid[j]+rfield.ConEmitReflec[j]) + 00458 rfield.reflin[j]); 00459 } 00460 } 00461 else if( nOption == 10 ) 00462 { 00463 /* this is exp(-tau) */ 00464 cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper ); 00465 /* need to free the vector once done */ 00466 lgFREE = true; 00467 for( j=0; j<nEnergy; ++j) 00468 { 00469 /* This is the TOTAL attenuation in both the continuum and lines. 00470 * Jon Miller discovered that the line attenuation was missing in 07.02 */ 00471 cont[j] = opac.ExpmTau[j]*rfield.trans_coef_total[j]; 00472 } 00473 } 00474 else 00475 { 00476 fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption); 00477 cdEXIT(EXIT_FAILURE); 00478 } 00479 00480 /* now generate the continua */ 00481 for( ncell = 0; ncell < nEnergy; ++ncell ) 00482 { 00483 if( ncell >= rfield.nflux ) 00484 { 00485 ReturnedSpectrum[ncell] = SMALLFLOAT; 00486 } 00487 else 00488 { 00489 ReturnedSpectrum[ncell] = cont[ncell]; 00490 } 00491 ASSERT( ReturnedSpectrum[ncell] >=0.f ); 00492 } 00493 00494 /* need to free the vector once done if this flag was set */ 00495 if( lgFREE ) 00496 { 00497 free(cont); 00498 } 00499 return; 00500 }