PLUTO
Functions | Variables
radiat.c File Reference
#include "pluto.h"
#include "cooling_defs.h"
Include dependency graph for radiat.c:

Go to the source code of this file.

Functions

void Radiat (double *v, double *rhs)
 
void Find_Rates (double T, double Ne, double N, double *v)
 
double find_N_rho ()
 
double H_MassFrac (void)
 

Variables

const double elem_ab [] = { 0.93, 0.074, 3.0e-4, 5.e-5, 4.0e-4, 7.0e-5, 1.5e-5, 2.69e-5 }
 
double elem_ab_est [] = { 0.93, 0.074, 3.0e-4, 5.e-5, 4.0e-4, 7.0e-5, 1.5e-5 }
 
double elem_ab_lod [] = { 0.92, 0.08, 2.3e-4, 6.2e-5, 4.4e-4, 7.0e-5, 1.25e-5 }
 
double elem_ab_sol [] = { 0.93, 0.074, 3.0e-4, 9.e-5, 7.e-4, 7.e-5, 1.e-5 }
 
double elem_ab_uni [] = { 0.93, 0.072, 5.0e-4, 9.e-5, 8.e-4, 8.e-5, 2.e-5 }
 
const double elem_mass [] = { 1.007, 4.002, 12.01, 14.01, 15.99, 20.18, 32.07, 55.845 }
 
const int elem_part []
 
const double rad_rec_z []
 
const double coll_ion_dE []
 
COOL_COEFF CoolCoeffs
 

Function Documentation

double find_N_rho ( )

Definition at line 694 of file radiat.c.

707 {
708  int i, j;
709  double mu, mu1, mu2;
710 
711  mu1 = mu2 = 0.0;
712  for (i = 0; i < 7; i++) {
713  mu1 += elem_ab[i]*elem_mass[i]; /* Numerator part of mu */
714  mu2 += elem_ab[i]; /* Denominator part of mu */
715  }
716  mu = mu1 / mu2;
717  return (UNIT_DENSITY/mu*CONST_NA); /* This is N/rho, with N the total number density of atoms and ions */
718 }
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
elem_mass
const double elem_mass[]
Definition: radiat.c:14
elem_ab
const double elem_ab[]
Definition: radiat.c:9
CONST_NA
#define CONST_NA
Avogadro Contant.
Definition: pluto.h:267
j
int j
Definition: analysis.c:2
i
int i
Definition: analysis.c:2

Here is the caller graph for this function:

void Find_Rates ( double  T,
double  Ne,
double  N,
double *  v 
)

Definition at line 413 of file radiat.c.

427 {
428  double dn, lam, t4, tmprec, ft1, ft2, tmpT, scrh;
429  int ti1, ti2, cindex, i, ions, nv, tindex, j, k;
430  static double ***ion_data, **intData;
431 
432  if (ion_data == NULL) { /* -- compute ionization rates tables -- */
433  ion_data = ARRAY_3D(I_g_stepNumber, 7, NIONS, double);
434  intData = ARRAY_2D(7, NIONS, double);
435  Create_Ion_Coeff_Tables(ion_data);
436  }
437 
438  for (nv = 0; nv < NIONS; nv++ ) {
439  CoolCoeffs.Rrate[nv] = 0.0;
440  CoolCoeffs.Lrate[nv] = 0.0;
441  CoolCoeffs.Crate[nv] = 0.0;
442  CoolCoeffs.Ra[nv] = 0.0;
443  CoolCoeffs.La[nv] = 0.0;
444  CoolCoeffs.Ca[nv] = 0.0;
445  CoolCoeffs.Rb[nv] = 0.0;
446  CoolCoeffs.Lb[nv] = 0.0;
447  CoolCoeffs.Cb[nv] = 0.0;
448  CoolCoeffs.Rc[nv] = 0.0;
449  CoolCoeffs.Lc[nv] = 0.0;
450  CoolCoeffs.Cc[nv] = 0.0;
451  }
452 
453  tmpT = T;
454  if (tmpT > I_TEND) tmpT = I_TEND - I_TSTEP*0.5;
455  else if (tmpT < I_TBEG) tmpT = I_TBEG + I_TSTEP*0.5;
456 
457  tindex = floor( (tmpT - I_TBEG)/I_TSTEP);
458 
459  ft1 = ( I_TBEG + I_TSTEP*((double)tindex + 1.0) - tmpT ) / I_TSTEP;
460  ft2 = ( tmpT - I_TBEG - I_TSTEP*((double)tindex) ) / I_TSTEP;
461 
462  for (k = 0; k < 7; k++) {
463  for (ions = 0; ions < NIONS; ions++) {
464  intData[k][ions] = ion_data[tindex][k][ions]*ft1 + ion_data[tindex + 1][k][ions]*ft2;
465  }}
466 
467  /* -------------------------------------------------------------
468  Computing rates for H
469  ------------------------------------------------------------- */
470 
471  lam = 157890. / T;
472  CoolCoeffs.Crate[0] += Ne*intData[0][el_H]; /* collisional ionization */
473  CoolCoeffs.Rrate[0] += Ne*intData[1][el_H]; /* total recombination - NIFS-DATA-54 - Kato & Asano 1999 , from Aldrovandi & Pequignot 1973 */
474 
475  CoolCoeffs.fCH = intData[0][el_H];
476  CoolCoeffs.fRH = intData[1][el_H];
477 
478  /* contribution from charge - exchange (recombination and ionization of all ions) : */
479 /*
480  for (ions=3; ions<28; ions++) {
481 
482  if (T/10000.<=chtrH_ion_upT[ions]) CoolCoeffs.Crate[0] += X[ions+1][ii] * N * elem_ab[elem_part[ions+1]] * 1.e-9 * chtrH_rec_a[ions] * pow( (T/10000.), chtrH_rec_b[ions] ) * ( 1. + chtrH_rec_c[ions] * exp ( chtrH_rec_d[ions]*(T/10000.) ));
483  if (T<=chtrH_rec_upT[ions]) CoolCoeffs.Rrate[0] += 1.e-9 * chtrH_ion_a[ions] * pow( (T/10000.), chtrH_ion_b[ions] ) * ( 1. + chtrH_ion_c[ions] * exp ( chtrH_ion_d[ions]*(T/10000.) ));
484 
485  tmprec = v[ions] * elem_ab[elem_part[ions]] / elem_ab[el_H] / fabs(v[HI]-1.e-13) * N * ion_data[4][ions];
486  CoolCoeffs.Crate[0] += tmprec;
487  CoolCoeffs.Ca[0] += ion_data[4][ions];
488  }
489 */
490 
491 /* -------------------------------------------------------------------
492  Do the other ions now
493  ------------------------------------------------------------------- */
494 
495  for (ions = 2; ions < NIONS - 1; ions++) {
496 
497  /* ------------------------------------------
498  collisional ionization - Voronov 1997
499  ------------------------------------------ */
500 
501  tmprec = Ne*intData[0][ions];
502  CoolCoeffs.Crate[ions] += tmprec;
503  CoolCoeffs.Lrate[ions+1] += tmprec;
504  CoolCoeffs.Ca[ions] += intData[0][ions];
505  CoolCoeffs.La[ions+1] += intData[0][ions];
506 
507  /* -------------------------------------------------------
508  radiative recombination - Pequignot & al, 1991 A&A
509  ------------------------------------------------------- */
510 
511  tmprec = Ne*intData[1][ions-1];
512  CoolCoeffs.Crate[ions] += tmprec;
513  CoolCoeffs.Rrate[ions-1] += tmprec;
514  CoolCoeffs.Ca[ions] += intData[1][ions-1];
515  CoolCoeffs.Ra[ions-1] += intData[1][ions-1];
516 
517  /* -----------------------------------------------------------
518  dielectronic recombination - Nussbaumer & Storey, 1983
519  ----------------------------------------------------------- */
520 
521  tmprec = Ne*intData[2][ions-1];
522  CoolCoeffs.Crate[ions] += tmprec;
523  CoolCoeffs.Rrate[ions-1] += tmprec;
524  CoolCoeffs.Ca[ions] += intData[2][ions-1];
525  CoolCoeffs.Ra[ions-1] += intData[2][ions-1];
526 
527  /* --------------------------------------------------------------------------
528  charge-transfer with H+ - ionization - Kingdon & Ferland 1996, ApJSS
529  -------------------------------------------------------------------------- */
530 
531  scrh = N*elem_ab[el_H]*intData[3][ions];
532  tmprec = (1.0 - v[X_HI])*scrh;
533  CoolCoeffs.Crate[ions] += tmprec;
534  CoolCoeffs.Lrate[ions+1] += tmprec;
535  CoolCoeffs.Cb[ions] -= scrh;
536  CoolCoeffs.Lb[ions+1] -= scrh;
537 
538  /* ---------------------------------------------------------------------------
539  charge-transfer with H - recombination - Kingdon & Ferland 1996, ApJSS
540  --------------------------------------------------------------------------- */
541 
542  scrh = N*elem_ab[el_H]*intData[4][ions-1];
543  tmprec = v[X_HI]*scrh;
544  CoolCoeffs.Crate[ions] += tmprec;
545  CoolCoeffs.Rrate[ions-1] += tmprec;
546  CoolCoeffs.Cb[ions] += scrh;
547  CoolCoeffs.Rb[ions-1] += scrh;
548  }
549 
550  /* -- calculation for last element (ions = NIONS - 1) -- */
551 
552  tmprec = Ne*intData[1][ions-1];
553  CoolCoeffs.Crate[ions] += tmprec;
554  CoolCoeffs.Rrate[ions-1] += tmprec;
555  CoolCoeffs.Ca[ions] += intData[1][ions-1];
556  CoolCoeffs.Ra[ions-1] += intData[1][ions-1];
557 
558  tmprec = Ne*intData[2][ions-1];
559  CoolCoeffs.Crate[ions] += tmprec;
560  CoolCoeffs.Rrate[ions-1] += tmprec;
561  CoolCoeffs.Ca[ions] += intData[2][ions-1];
562  CoolCoeffs.Ra[ions-1] += intData[2][ions-1];
563 
564  scrh = N*elem_ab[el_H]*intData[4][ions-1];
565  tmprec = v[X_HI]*scrh;
566  CoolCoeffs.Crate[ions] += tmprec;
567  CoolCoeffs.Rrate[ions-1] += tmprec;
568  CoolCoeffs.Cb[ions] += scrh;
569  CoolCoeffs.Rb[ions-1] += scrh;
570 
571  /* -- calculation for HeI -- */
572 
573  ions = 1;
574  tmprec = Ne*intData[0][ions];
575  CoolCoeffs.Crate[ions] += tmprec;
576  CoolCoeffs.Lrate[ions+1] += tmprec;
577  CoolCoeffs.Ca[ions] += intData[0][ions];
578  CoolCoeffs.La[ions+1] += intData[0][ions];
579 
580  scrh = N*elem_ab[el_H]*intData[3][ions];
581  tmprec = (1.0 - v[X_HI])*scrh;
582  CoolCoeffs.Crate[ions] += tmprec;
583  CoolCoeffs.Lrate[ions+1] += tmprec;
584  CoolCoeffs.Cb[ions] -= scrh;
585  CoolCoeffs.Lb[ions+1] -= scrh;
586 
587 /* --------------------------------------------------------------------------
588  charge-transfer with He - recombination
589  ORNL Charge Transfer Database
590  http://www-cfadc.phy.ornl.gov/astro/ps/data/cx/helium/rates/fits.data
591  -------------------------------------------------------------------------- */
592 
593  for (ions = 4; ions < NIONS; ions++) { /* 4 ??????????????? */
594  scrh = N*elem_ab[el_He]*intData[5][ions-1];
595  tmprec = v[X_HeI]*scrh;
596  CoolCoeffs.Crate[ions] += tmprec;
597  CoolCoeffs.Rrate[ions-1] += tmprec;
598 
599  CoolCoeffs.Cc[ions] += scrh;
600  CoolCoeffs.Rc[ions-1] += scrh;
601  }
602 
603 /* -------------------------------------------------------
604  And now, for ions for which we only have the total
605  electron-ion recombination coefficient:
606  ( NIFS-DATA-54 - Kato & Asano 1999 )
607  ------------------------------------------------------- */
608 
609  for (ions = 2; ions < NIONS; ions++) {
610  tmprec = Ne*intData[6][ions-1];
611  CoolCoeffs.Crate[ions] += tmprec;
612  CoolCoeffs.Rrate[ions-1] += tmprec;
613 
614  /* -----------------------------------------------
615  Now save the data for the partial derivatives
616  computation
617  ----------------------------------------------- */
618 
619  CoolCoeffs.Ca[ions] += intData[6][ions-1];
620  CoolCoeffs.Ra[ions-1] += intData[6][ions-1];
621  }
622 
623  /* He/He+ contribution from charge-exchange processes: */
624 /*
625  for (ions=3; ions<28; ions++) {
626  if (T<5000.) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a1[ions] * pow( (T/10000.), chtrHe_rec_b1[ions] ) * ( 1. + chtrHe_rec_c1[ions] * exp ( chtrHe_rec_d1[ions]*(T/10000.) ));
627  if ( (T>=5000.) && (T<10000.) ) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a2[ions] * pow( (T/10000.), chtrHe_rec_b2[ions] ) * ( 1. + chtrHe_rec_c2[ions] * exp ( chtrHe_rec_d2[ions]*(T/10000.) ));
628  if ( (T>=10000.) && (T<50000.) ) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a3[ions] * pow( (T/10000.), chtrHe_rec_b3[ions] ) * ( 1. + chtrHe_rec_c3[ions] * exp ( chtrHe_rec_d3[ions]*(T/10000.) ));
629  if ( (T>=50000.) && (T<100000.) ) CoolCoeffs.Rrate[2] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a4[ions] * pow( (T/10000.), chtrHe_rec_b4[ions] ) * ( 1. + chtrHe_rec_c4[ions] * exp ( chtrHe_rec_d4[ions]*(T/10000.) ));
630  if (T>=100000.) CoolCoeffs.Rrate[2] += N * v[ions] *elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a5[ions] * pow( (T/10000.), chtrHe_rec_b5[ions] ) * ( 1. + chtrHe_rec_c5[ions] * exp ( chtrHe_rec_d5[ions]*(T/10000.) ));
631 
632  if (T<5000.) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a1[ions] * pow( (T/10000.), chtrHe_rec_b1[ions] ) * ( 1. + chtrHe_rec_c1[ions] * exp ( chtrHe_rec_d1[ions]*(T/10000.) ));
633  if ( (T>=5000.) && (T<10000.) ) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a2[ions] * pow( (T/10000.), chtrHe_rec_b2[ions] ) * ( 1. + chtrHe_rec_c2[ions] * exp ( chtrHe_rec_d2[ions]*(T/10000.) ));
634  if ( (T>=10000.) && (T<50000.) ) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a3[ions] * pow( (T/10000.), chtrHe_rec_b3[ions] ) * ( 1. + chtrHe_rec_c3[ions] * exp ( chtrHe_rec_d3[ions]*(T/10000.) ));
635  if ( (T>=50000.) && (T<100000.) ) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a4[ions] * pow( (T/10000.), chtrHe_rec_b4[ions] ) * ( 1. + chtrHe_rec_c4[ions] * exp ( chtrHe_rec_d4[ions]*(T/10000.) ));
636  if (T>=100000.) CoolCoeffs.Crate[1] += N * v[ions] * elem_ab[elem_part[ions]] * 1.e-9 * chtrHe_rec_a5[ions] * pow( (T/10000.), chtrHe_rec_b5[ions] ) * ( 1. + chtrHe_rec_c5[ions] * exp ( chtrHe_rec_d5[ions]*(T/10000.) ));
637  }
638 */
639 
640 /* -------------------------------------------------------------
641  C H E C K
642  ------------------------------------------------------------- */
643 /*
644 for (ions = 0; ions < NIONS; ions++) {
645  CoolCoeffs.Lrate[ions] *= UNIT_LENGTH/UNIT_VELOCITY;
646  CoolCoeffs.Crate[ions] *= UNIT_LENGTH/UNIT_VELOCITY;
647  CoolCoeffs.Rrate[ions] *= UNIT_LENGTH/UNIT_VELOCITY;
648 }
649 
650 print ("H: %12.6e %12.6e %12.6e\n",
651  CoolCoeffs.Lrate[HI-NFLX], CoolCoeffs.Crate[HI-NFLX],CoolCoeffs.Rrate[HI-NFLX]);
652 
653 print ("--------------------------------------------------------------- \n");
654 for (nv = HeI; nv <= HeII; nv++){
655 print ("He(%d): %12.6e %12.6e %12.6e\n",nv-HeI+1,
656  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
657 }
658 
659 print ("--------------------------------------------------------------- \n");
660 for (nv = CI; nv <= CV; nv++){
661 print ("C(%d): %12.6e %12.6e %12.6e\n",nv-CI+1,
662  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
663 }
664 print ("--------------------------------------------------------------- \n");
665 
666 for (nv = NI; nv <= NV; nv++){
667 print ("N(%d): %12.6e %12.6e %12.6e\n",nv-NI+1,
668  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
669 }
670 print ("--------------------------------------------------------------- \n");
671 
672 for (nv = OI; nv <= OV; nv++){
673 print ("O(%d): %12.6e %12.6e %12.6e\n",nv-OI+1,
674  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
675 }
676 print ("--------------------------------------------------------------- \n");
677 
678 for (nv = NeI; nv <= NeV; nv++){
679 print ("Ne(%d): %12.6e %12.6e %12.6e\n",nv-NeI+1,
680  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
681 }
682 print ("--------------------------------------------------------------- \n");
683 for (nv = SI; nv <= SV; nv++){
684 print ("Se(%d): %12.6e %12.6e %12.6e\n",nv-SI+1,
685  CoolCoeffs.Lrate[nv-NFLX], CoolCoeffs.Crate[nv-NFLX],CoolCoeffs.Rrate[nv-NFLX]);
686 }
687 
688 
689 exit(1);
690 */
691 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
I_TBEG
#define I_TBEG
Definition: cooling_defs.h:171
COOLING_COEFFICIENTS::Crate
double Crate[NIONS]
Definition: cooling_defs.h:31
COOLING_COEFFICIENTS::Cc
double Cc[NIONS]
Definition: cooling_defs.h:41
COOLING_COEFFICIENTS::Rb
double Rb[NIONS]
Definition: cooling_defs.h:43
configure.scrh
tuple scrh
Definition: configure.py:200
COOLING_COEFFICIENTS::Lrate
double Lrate[NIONS]
Definition: cooling_defs.h:30
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
COOLING_COEFFICIENTS::Lb
double Lb[NIONS]
Definition: cooling_defs.h:37
el_He
#define el_He
Definition: cooling_defs.h:72
COOLING_COEFFICIENTS::fRH
double fRH
Definition: cooling_defs.h:52
elem_ab
const double elem_ab[]
Definition: radiat.c:9
COOLING_COEFFICIENTS::La
double La[NIONS]
Definition: cooling_defs.h:36
I_TSTEP
#define I_TSTEP
Definition: cooling_defs.h:170
X_HeI
Definition: cooling.h:110
I_TEND
#define I_TEND
Definition: cooling_defs.h:172
NIONS
#define NIONS
Definition: cooling.h:28
I_g_stepNumber
#define I_g_stepNumber
Definition: cooling_defs.h:169
j
int j
Definition: analysis.c:2
X_HI
Definition: cooling.h:110
k
int k
Definition: analysis.c:2
COOLING_COEFFICIENTS::fCH
double fCH
Definition: cooling_defs.h:52
COOLING_COEFFICIENTS::Rrate
double Rrate[NIONS]
Definition: cooling_defs.h:32
COOLING_COEFFICIENTS::Rc
double Rc[NIONS]
Definition: cooling_defs.h:44
i
int i
Definition: analysis.c:2
COOLING_COEFFICIENTS::Ra
double Ra[NIONS]
Definition: cooling_defs.h:42
COOLING_COEFFICIENTS::Lc
double Lc[NIONS]
Definition: cooling_defs.h:38
COOLING_COEFFICIENTS::Ca
double Ca[NIONS]
Definition: cooling_defs.h:39
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
Create_Ion_Coeff_Tables
int Create_Ion_Coeff_Tables(double ***)
Definition: make_tables.c:200
COOLING_COEFFICIENTS::Cb
double Cb[NIONS]
Definition: cooling_defs.h:40
CoolCoeffs
COOL_COEFF CoolCoeffs
Definition: radiat.c:38
el_H
#define el_H
Definition: cooling_defs.h:71

Here is the call graph for this function:

Here is the caller graph for this function:

double H_MassFrac ( void  )

Compute the mass fraction X of Hydrogen as function of the composition of the gas. 

        f_H A_H

X = -------------— f_k A_k

where

f_k : is the fractional abundance (by number), f_k = N_k/N_tot of atomic species (no matter ionization degree).

A_K : is the atomic weight

where N_{tot} is the total number density of particles

ARGUMENTS

none

Definition at line 721 of file radiat.c.

746 {
747  double mmw2;
748  int i, j;
749 
750  mmw2 = 0.0;
751  for (i = 0; i < (Fe_IONS==0?7:8); i++) {
752  mmw2 += elem_mass[i]*elem_ab[i]; /* Denominator part of X */
753  }
754  return ( (elem_mass[0]*elem_ab[0]) / mmw2);
755 }
elem_mass
const double elem_mass[]
Definition: radiat.c:14
elem_ab
const double elem_ab[]
Definition: radiat.c:9
j
int j
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
Fe_IONS
#define Fe_IONS
Definition: cooling.h:18
void Radiat ( double *  v,
double *  rhs 
)

Cooling for optically thin plasma up to about 200,000 K Plasma composition: H, HeI-II, CI-V, NI-V, OI-V, NeI-V, SI-V Assumed abundances in elem_ab Uses S : Array = Variables vector x line points rhs : output for the system of ODE ibeg, iend : begin and end points of the current line

Definition at line 40 of file radiat.c.

51 {
52  int nv, j, k, cooling_nT, cooling_nNe, jj;
53  int ti1, ti2, ne1, ne2, nrt;
54  double mu, sT, scrh, tmpT, tmpNe, tt1, tt2, nn1, nn2, tf1, tf2, nf1, nf2;
55  double N, n_el, rlosst, em, cf1, cf2;
56  double T, *X, *RS;
57  static double ***tab;
58  static double N_rho;
59  static double E_cost, Unit_Time; /* -- for dimensionalization purposes -- */
60  double em2, em3;
61 double prs;
62 
63 /* --------------------------------------------------------------------
64  Load tables from disk at first call.
65  Tables are 2-D with T and Ne being the coordinates.
66  -------------------------------------------------------------------- */
67 
68  if (tab == NULL) {
69 
70  E_cost = UNIT_LENGTH/UNIT_DENSITY/
71  (UNIT_VELOCITY*UNIT_VELOCITY*UNIT_VELOCITY);
72  Unit_Time = UNIT_LENGTH/UNIT_VELOCITY;
73 
74  /* ---------------------------------------
75  Compute cooling function tables
76  --------------------------------------- */
77 
78  for (nv = 0; nv < NIONS; nv++) CoolCoeffs.dLIR_dX[nv] = 0.0;
79 
80  n_el = C_NeMIN;
81  ne1 = 0;
82  while (n_el < C_NeMAX) {
83  T = C_TMIN;
84  ne2 = 0;
85  while (T < C_TMAX) {
86  tmpT = T;
87  T = T*exp(C_TSTEP); /* should be *exp(0.02) */
88  ne2 = ne2 + 1;
89  }
90  tmpNe = n_el;
91  n_el = n_el*exp(C_NeSTEP); /* should be *exp(0.06) */
92  ne1 = ne1 + 1;
93  }
94 
95  tab = ARRAY_3D(NIONS, ne1, ne2, double);
96  Create_Losses_Tables(tab, &cooling_nT, &cooling_nNe);
97  N_rho = find_N_rho();
98  } /* end load tables */
99 
100 /* ---------------------------------------
101  Force species to lie between 0 and 1
102  --------------------------------------- */
103 
104  for (nv = NFLX; nv < (NFLX + NIONS); nv++){
105  v[nv] = MIN(v[nv], 1.0);
106  v[nv] = MAX(v[nv], 0.0);
107  }
108 
109  prs = v[RHOE]*(g_gamma-1.0);
110  if (prs < 0.0) {
111  prs = g_smallPressure;
112  v[RHOE] = prs/(g_gamma - 1.0);
113  }
114 
115  mu = MeanMolecularWeight(v);
116  T = prs/v[RHO]*KELVIN*mu;
117 
118  if (mu < 0.0){
119  print ("! Radiat: negative mu\n");
120  QUIT_PLUTO(1);
121  }
122 
123  /* ---------------------------------
124  offset pointers to work more
125  efficiently.
126  --------------------------------- */
127 
128  X = v + NFLX;
129  RS = rhs + NFLX;
130 
131  sT = sqrt(T);
132  N = v[RHO]*N_rho; /* -- Total number density -- */
133 
134 /* -----------------------------------
135  compute electron number density
136  ----------------------------------- */
137 
138  CoolCoeffs.dnel_dX[0] = -N*elem_ab[el_H];
139  n_el = N*(1.0 - v[X_HI])*elem_ab[el_H]; /* -- contribution from ionized H -- */
140 
141  for (nv = 1; nv < NIONS; nv++) {
142  CoolCoeffs.dnel_dX[nv] = N*(rad_rec_z[nv] - 1.0)*elem_ab[elem_part[nv]];
143  n_el += X[nv]*CoolCoeffs.dnel_dX[nv];
144  }
145 
146  if (n_el/N < 1.e-4) n_el = N*1.e-4; /* OK ????? */
147 
148  CoolCoeffs.Ne = n_el;
149 
150  Find_Rates(T, n_el, N, v); /* find transition rates */
151 
152 /* -----------------------------------------------------------------------------
153  Compute RHS of the system of ODE.
154  Evaluate the transition rates
155  (multiply the rates given by find_rates()
156  by numerical densities)
157  ----------------------------------------------------------------------------- */
158 
159  /* ----------------------------------------------
160  Compute right hand side for all ions
161  ---------------------------------------------- */
162 
163  for (nv = 1; nv <= NIONS - 2; nv++) {
164  RS[nv] = X[nv - 1]*CoolCoeffs.Lrate[nv]
165  - X[nv] *CoolCoeffs.Crate[nv]
166  + X[nv + 1]*CoolCoeffs.Rrate[nv];
167  }
168  RS[0] = (1.0 - X[0])*CoolCoeffs.Rrate[0] - X[0]*CoolCoeffs.Crate[0]; /* separate for H */
169 
170  nv = NIONS - 1;
171  RS[nv] = X[nv - 1]*CoolCoeffs.Lrate[nv] - X[nv]*CoolCoeffs.Crate[nv];
172 
173  for (nv = 0; nv < NIONS; nv++) RS[nv] *= Unit_Time;
174 
175 /* ---------------------------------------------------------------
176  check sums
177  --------------------------------------------------------------- */
178 /*
179 {
180 double sum;
181 
182 sum=fabs(RS[1] + RS[2]);
183 if (sum > 1.e-17){
184  print("> Sum(RHS) != 0 for He!!\n");
185  exit(1);
186 }
187 
188 sum = 0.0;
189 for (nv = 0; nv < C_IONS; nv++) sum += RS[CI-NFLX+nv];
190 if (fabs(sum) > 1.e-10){
191  print("> Sum(RHS) != 0 for C !! (%12.6e)\n", sum);
192  exit(1);
193 }
194 }
195 
196 nv = 8;
197 sum=fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
198 if ( sum > 1.e-10){
199  print("> Sum(RHS) != 0 for N!!\n");
200  exit(1);
201 }
202 nv = 13;
203 sum=fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
204 if ( sum > 1.e-10){
205  print("> Sum(RHS) != 0 for O!!\n");
206  exit(1);
207 }
208 nv = 18;
209 sum = fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
210 if ( sum > 1.e-10){
211  print("> Sum(RHS) != 0 for Ne!!\n");
212  exit(1);
213 }
214 nv = 23;
215 fabs(RS[nv] + RS[nv+1] + RS[nv+2] + RS[nv+3] + RS[nv+4]);
216 if ( sum > 1.e-10){
217  print("> Sum(RHS) != 0 for S!!\n");
218  exit(1);
219 }
220 }
221 */
222 
223 /*
224 print("%12.6e %12.6e %12.6e %12.6e %12.6e %12.6e\n", v[HI],
225  RS[0],RS[1], RS[2], RS[3], RS[4]);
226 }
227 exit(1);
228 */
229 
230 /* ********************************************************************************************************* */
231 
232 
233  /* -------------------------------------------------------------
234  Find the indexes in the tables needed for interpolation
235 
236  Constrain both T and n_el to lie in the range provided
237  by the tables.
238  This does not affect the original vector of primitive
239  quantities, but it is more computationally efficient.
240  In other words, when T > Tmax, the cooling function
241  will saturate at Tmax.
242  ------------------------------------------------------------- */
243 
244  tmpT = T; /* keep the double value for after the interpolation */
245  tmpNe = n_el;
246 
247  if (tmpT <= C_TMIN) tmpT = C_TMIN + C_TSTEP*0.5;
248  if (tmpT >= C_TMAX) tmpT = C_TMAX - C_TSTEP*0.5;
249 
250  if (n_el <= C_NeMIN) tmpNe = C_NeMIN + C_NeSTEP*0.5;
251  if (n_el >= C_NeMAX) tmpNe = C_NeMAX - C_NeSTEP*0.5;
252 
253  scrh = log(tmpT/C_TMIN)/C_TSTEP;
254  ti1 = floor(scrh);
255  ti2 = ceil (scrh);
256 
257  scrh = log(tmpNe/C_NeMIN)/C_NeSTEP;
258  ne1 = floor(scrh);
259  ne2 = ceil (scrh);
260 
261  tt1 = C_TMIN*exp(ti1*C_TSTEP);
262  tt2 = C_TMIN*exp(ti2*C_TSTEP);
263 
264  nn1 = C_NeMIN*exp(ne1*C_NeSTEP);
265  nn2 = C_NeMIN*exp(ne2*C_NeSTEP);
266 
267  tf1 = (tt2 - tmpT)/(tt2 - tt1);
268  tf2 = (tmpT - tt1)/(tt2 - tt1);
269  nf1 = (nn2 - tmpNe)/(nn2 - nn1);
270  nf2 = (tmpNe - nn1)/(nn2 - nn1);
271 
272 /* ------------------------------------------------------------
273  get \Delta e[k] coefficients
274  ------------------------------------------------------------ */
275 
276  for (nv = 0; nv < NIONS; nv++) {
277  cf1 = tf1*tab[nv][ne1][ti1] + tf2*tab[nv][ne1][ti2];
278  cf2 = tf1*tab[nv][ne2][ti1] + tf2*tab[nv][ne2][ti2];
279  CoolCoeffs.de[nv] = nf1*cf1 + nf2*cf2;
280  }
281 
282 /* -- Save the slopes of the radiative losses coefficients from the table -- */
283 
284  scrh = 0.5/(nn2 - nn1);
285  for (nv = 0; nv < NIONS; nv++) {
286  CoolCoeffs.de_dne[nv] = ( tab[nv][ne2][ti1] + tab[nv][ne2][ti2]
287  - tab[nv][ne1][ti1] - tab[nv][ne1][ti2])*scrh;
288  }
289 
290 /* ---------------------------------------------------------------
291  For HeII we impose a temperature cut-off above 9.e4 K
292  since HeII --> HeIII
293  --------------------------------------------------------------- */
294 
295  if (T >= 9.e4) {
296  scrh = exp(-(T - 9.e+4)/4.e+4);
297  CoolCoeffs.de[2] *= scrh;
298  CoolCoeffs.de_dne[2] *= scrh;
299  }
300 
301  /* ... and then continue with the radiative losses computation */
302 
303  em = 0.0;
304  for (nv = 0; nv < NIONS; nv++) {
305  em += X[nv]*CoolCoeffs.de[nv]*elem_ab[elem_part[nv]];
306  }
307 
308  if (em < 0.0) print("> 2 negative em => possitive losses???? em = %12.6e\n",em);
309 
310 /* -------------------------------------------------------------------
311  Now finally add the ionization / recombination losses ...
312  ------------------------------------------------------------------- */
313 
314  /* -----------------------------------
315  thermal energy lost by ionization
316  ----------------------------------- */
317 
318  CoolCoeffs.dLIR_dX[0] = 1.08e-8*sT*elem_ab[el_H]/13.6*exp(-157890.0/T)*CONST_eV;
319  em += CoolCoeffs.dLIR_dX[0]*v[X_HI];
320 
321  /*
322  if ( (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV) / (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV) > 2.0) {
323  print("too different em's: %12.6e - %12.6e\n",
324  (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV), (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV));
325 
326  }
327  if ( (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV) / (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV) < 0.5) {
328  print("too different em's: %12.6e - %12.6e\n",
329  (1.08e-8*sT*v[HI]*elem_ab[el_H]/13.6*exp(-157890.0/T) *CONST_eV), (v[HI] *CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*Unit_Time*CONST_eV));
330 
331  }
332 */
333 /*
334  em2 = v[HI]*CoolCoeffs.Crate[0]*elem_ab[el_H]*coll_ion_dE[0]*CONST_eV*N;
335  for (nv=2; nv<NIONS; nv++) {
336  if (CoolCoeffs.Lrate[nv]>0) em2 += X[nv-1]*CoolCoeffs.Lrate[nv]*elem_ab[elem_part[nv]]*coll_ion_dE[nv]*CONST_eV*N;
337  }
338 */
339  /* --------------------------------------
340  thermal energy lost by recombination
341  (2/3 kT per recombination).
342  -------------------------------------- */
343 
344  scrh = 2.6e-11/sT*elem_ab[el_H]*T/11590.0*0.67*CONST_eV;
345  em += scrh*(1.0 - v[X_HI]);
346  CoolCoeffs.dLIR_dX[0] -= scrh;
347 
348 /*
349  if (T>4000) print("ratesHI: old %12.6e new %12.6e at T = %12.6e\n",2.6e-11/sT*(1.0 - v[HI])*elem_ab[el_H]*T/11590.0*0.67*CONST_eV,
350  (1.0 - v[HI])*CoolCoeffs.Rrate[0]*elem_ab[el_H]*0.67*CONST_kB*T,T);
351  */
352 /*
353  em3 = (1.0 - v[HI])*CoolCoeffs.Rrate[0]*elem_ab[el_H]*0.67*CONST_kB*T*N;
354  for (nv=1; nv<NIONS-1; nv++) {
355  if (CoolCoeffs.Rrate[nv]>0) em3 += X[nv+1]*CoolCoeffs.Rrate[nv]*elem_ab[elem_part[nv]]*0.67*CONST_kB*T*N;
356  }
357 */
358 
359  /* --------------------------------------------------------
360  Add cooling contribution from bremsstrahlung
361  -------------------------------------------------------- */
362 
363  scrh = 1.42e-27*sT*elem_ab[el_H];
364  CoolCoeffs.dLIR_dX[2] = 1.42e-27*sT*elem_ab[el_He];
365  em += scrh*(1.0 - v[X_HI]) + CoolCoeffs.dLIR_dX[2]*v[X_HeII];
366 
367  CoolCoeffs.dLIR_dX[0] -= scrh;
368 
369  /* --------------------------------------------------------------------
370  Add cooling contribution from FeII 1.6 & 25 micron
371  and Mg II 2800 and Si II 35 micron
372  emission lines. This is TEMPORARY until the adding
373  of FeI and FeII as ions to NEq !
374  -------------------------------------------------------------------- */
375 
376  em += 1.6e-12*8.63e-6*0.3*0.0495*exp(-575.0/T)/sT * 580.0*sT/(n_el + 580.0*sT) *0.00004;
377  em += 1.6e-12*8.63e-6*0.39*0.775*exp(-8980.0/T)/sT * 1130.0*sT/(n_el + 1130.0*sT) *0.00004;
378 
379  em += 1.6e-12*8.63e-6*8.0*4.43*exp(-5.13e4/T)/sT*1.e10*sT/(n_el + 1.e10*sT)*0.00002*0.7; /* Mg II 2800: Mendoza */
380  em += 1.6e-12*8.63e-6*2.85*0.0354*exp(-410.0/T)/sT*16.8*sT/(n_el + 16.8*sT)*0.000033*0.7; /* Si II 35 micron: Dufton&Kingston, MNRAS 248 */
381 
382  /* ---------------------------------------------------
383  rlosst is the energy loss in units of erg/cm^3/s;
384  it must be multiplied by cost_E in order to match
385  non-dimensional units.
386  Source term for the neutral fraction scales with
387  Unit_Time
388 
389  Pressure source term equals zero when
390  T < g_minCoolingTemp
391  --------------------------------------------------- */
392 
393  rlosst = em*n_el*N;
394 /*
395  rlosst += em3 + em2;
396 */
397 
398 /*
399  if (T > g_minCoolingTemp) rhs[PRS] = -E_cost*rlosst*(g_gamma - 1.0);
400  else rhs[PRS] = 0.0;
401 */
402 
403  rhs[RHOE] = -E_cost*rlosst;
404  rhs[RHOE] *= 1.0/(1.0 + exp( -(T - g_minCoolingTemp)/100.)); /* -- cut cooling function --*/
405 
406  if (rhs[RHOE] > 0.0) {
407  print ("! Error: positive radiative losses! %12.6e %12.6e %12.6e \n", em, n_el, N);
408  QUIT_PLUTO(1);
409  }
410 }
COOLING_COEFFICIENTS::dnel_dX
double dnel_dX[NIONS]
Definition: cooling_defs.h:51
MAX
#define MAX(a, b)
Definition: macros.h:101
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
COOLING_COEFFICIENTS::dLIR_dX
double dLIR_dX[NIONS]
Definition: cooling_defs.h:55
COOLING_COEFFICIENTS::Crate
double Crate[NIONS]
Definition: cooling_defs.h:31
g_gamma
double g_gamma
Definition: globals.h:112
COOLING_COEFFICIENTS::de_dne
double de_dne[NIONS]
Definition: cooling_defs.h:47
UNIT_DENSITY
#define UNIT_DENSITY
Unit density in gr/cm^3.
Definition: pluto.h:369
CONST_eV
#define CONST_eV
Electron Volt in erg.
Definition: pluto.h:256
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
COOLING_COEFFICIENTS::de
double de[NIONS]
Definition: cooling_defs.h:46
COOLING_COEFFICIENTS::Lrate
double Lrate[NIONS]
Definition: cooling_defs.h:30
C_TSTEP
#define C_TSTEP
Definition: cooling_defs.h:185
g_smallPressure
double g_smallPressure
Small value for pressure fix.
Definition: globals.h:110
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
KELVIN
#define KELVIN
Definition: pluto.h:401
el_He
#define el_He
Definition: cooling_defs.h:72
elem_ab
const double elem_ab[]
Definition: radiat.c:9
MIN
#define MIN(a, b)
Definition: macros.h:104
rad_rec_z
const double rad_rec_z[]
Definition: radiat.c:24
UNIT_VELOCITY
#define UNIT_VELOCITY
Unit velocity in cm/sec.
Definition: pluto.h:377
elem_part
const int elem_part[]
Definition: radiat.c:16
C_TMAX
#define C_TMAX
Definition: cooling_defs.h:184
find_N_rho
double find_N_rho()
Definition: radiat.c:694
X
#define X
Definition: radiat.modified_by_Bhargav.c:3
NFLX
#define NFLX
Definition: mod_defs.h:32
UNIT_LENGTH
#define UNIT_LENGTH
Unit Length in cm.
Definition: pluto.h:373
NIONS
#define NIONS
Definition: cooling.h:28
C_NeSTEP
#define C_NeSTEP
Definition: cooling_defs.h:181
j
int j
Definition: analysis.c:2
X_HI
Definition: cooling.h:110
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
COOLING_COEFFICIENTS::Rrate
double Rrate[NIONS]
Definition: cooling_defs.h:32
g_minCoolingTemp
double g_minCoolingTemp
The minimum temperature (in K) below which cooling is suppressed.
Definition: globals.h:106
C_TMIN
#define C_TMIN
Definition: cooling_defs.h:183
MeanMolecularWeight
double MeanMolecularWeight(double *V)
Definition: radiat.modified_by_Bhargav.c:178
Find_Rates
void Find_Rates(double T, double Ne, double N, double *v)
Definition: radiat.c:413
C_NeMAX
#define C_NeMAX
Definition: cooling_defs.h:180
Create_Losses_Tables
int Create_Losses_Tables(double ***, int *, int *)
Definition: make_tables.c:841
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
CoolCoeffs
COOL_COEFF CoolCoeffs
Definition: radiat.c:38
el_H
#define el_H
Definition: cooling_defs.h:71
COOLING_COEFFICIENTS::Ne
double Ne
Definition: cooling_defs.h:34
C_NeMIN
#define C_NeMIN
Definition: cooling_defs.h:179

Here is the call graph for this function:

Variable Documentation

const double coll_ion_dE[]
Initial value:
= { 13.6, 24.6, 54.4
Ne_EXPAND(21.6, 41.0, 63.5, 97.1, 126.2)
Fe_EXPAND(7.87, 16.1879, 30.652) }
Ne_EXPAND
Definition: cooling.h:114
Fe_EXPAND
#define Fe_EXPAND(a, b, c)
Definition: cooling.h:91

Definition at line 31 of file radiat.c.

COOL_COEFF CoolCoeffs

Definition at line 38 of file radiat.c.

const double elem_ab[] = { 0.93, 0.074, 3.0e-4, 5.e-5, 4.0e-4, 7.0e-5, 1.5e-5, 2.69e-5 }

Definition at line 9 of file radiat.c.

double elem_ab_est[] = { 0.93, 0.074, 3.0e-4, 5.e-5, 4.0e-4, 7.0e-5, 1.5e-5 }

Definition at line 10 of file radiat.c.

double elem_ab_lod[] = { 0.92, 0.08, 2.3e-4, 6.2e-5, 4.4e-4, 7.0e-5, 1.25e-5 }

Definition at line 11 of file radiat.c.

double elem_ab_sol[] = { 0.93, 0.074, 3.0e-4, 9.e-5, 7.e-4, 7.e-5, 1.e-5 }

Definition at line 12 of file radiat.c.

double elem_ab_uni[] = { 0.93, 0.072, 5.0e-4, 9.e-5, 8.e-4, 8.e-5, 2.e-5 }

Definition at line 13 of file radiat.c.

const double elem_mass[] = { 1.007, 4.002, 12.01, 14.01, 15.99, 20.18, 32.07, 55.845 }

Definition at line 14 of file radiat.c.

const int elem_part[]
Initial value:
= {0, 1, 1
Ne_EXPAND(5, 5, 5, 5, 5)
Fe_EXPAND(7, 7, 7) }
Ne_EXPAND
Definition: cooling.h:114
Fe_EXPAND
#define Fe_EXPAND(a, b, c)
Definition: cooling.h:91

Definition at line 16 of file radiat.c.

const double rad_rec_z[]
Initial value:
= { 1., 1., 2.
Ne_EXPAND(1., 2., 3., 4., 5.)
Fe_EXPAND(1., 2., 3.) }
Ne_EXPAND
Definition: cooling.h:114
Fe_EXPAND
#define Fe_EXPAND(a, b, c)
Definition: cooling.h:91

Definition at line 24 of file radiat.c.