00001
00002
00003
00004 #include "cddefines.h"
00005 #include "ipoint.h"
00006 #include "rfield.h"
00007 #include "two_photon.h"
00008
00009 void TwoPhotonSetup( vector<two_photon> &tnu_vec, const long &ipHi, const long &ipLo, const double &Aul, const TransitionProxy &tr, const long ipISO, const long nelem )
00010 {
00011 DEBUG_ENTRY( "TwoPhotonSetup()" );
00012
00013 tnu_vec.resize( tnu_vec.size() + 1 );
00014 two_photon &tnu = tnu_vec.back();
00015
00016 tnu.ipHi = ipHi;
00017 tnu.ipLo = ipLo;
00018 tnu.AulTotal = Aul;
00019 tnu.Pop = &(*tr.Hi()).Pop();
00020 tnu.E2nu = tr.EnergyRyd();
00021 tnu.induc_dn_max = 0.;
00022
00023
00024 tnu.ipTwoPhoE = ipoint(tnu.E2nu);
00025 while( rfield.anu[tnu.ipTwoPhoE] > tnu.E2nu )
00026 {
00027 --tnu.ipTwoPhoE;
00028 }
00029 tnu.ipSym2nu.resize( tnu.ipTwoPhoE );
00030 tnu.As2nu.resize( tnu.ipTwoPhoE );
00031 tnu.local_emis.resize( tnu.ipTwoPhoE );
00032
00033
00034 for( long i=0; i < tnu.ipTwoPhoE; i++ )
00035 {
00036
00037 double energy = tnu.E2nu - rfield.anu[i];
00038
00039
00040 energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
00041
00042 tnu.ipSym2nu[i] = ipoint(energy);
00043 while( rfield.anu[tnu.ipSym2nu[i]] > tnu.E2nu ||
00044 tnu.ipSym2nu[i] >= tnu.ipTwoPhoE)
00045 {
00046 --tnu.ipSym2nu[i];
00047 }
00048 ASSERT( tnu.ipSym2nu[i] >= 0 );
00049 }
00050
00051 double SumShapeFunction = 0., Renorm= 0.;
00052
00053
00054
00055 for( long i=0; i < tnu.ipTwoPhoE; i++ )
00056 {
00057 double ShapeFunction;
00058
00059 ASSERT( rfield.anu[i]<=tnu.E2nu );
00060 ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/tnu.E2nu, ipISO, nelem )*rfield.widflx[i]/tnu.E2nu;
00061 SumShapeFunction += ShapeFunction;
00062
00063
00064
00065 tnu.As2nu[i] = (realnum)( tnu.AulTotal * ShapeFunction );
00066 }
00067
00068
00069
00070 Renorm = 1./SumShapeFunction;
00071
00072 for( long i=0; i < tnu.ipTwoPhoE; i++ )
00073 {
00074 tnu.As2nu[i] *= (realnum)Renorm;
00075 }
00076
00077
00078 ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
00079
00080 return;
00081 }
00082
00083 void CalcTwoPhotonRates( two_photon& tnu, bool lgDoInduced )
00084 {
00085 DEBUG_ENTRY( "CalcTwoPhotonRates()" );
00086
00087
00088 ASSERT( tnu.ipTwoPhoE>0 && tnu.E2nu>0. );
00089
00090 double sum = 0.;
00091 tnu.induc_up = 0.;
00092 tnu.induc_dn = 0.;
00093
00094
00095 for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
00096 {
00097
00098
00099 ASSERT( rfield.anu[nu] < 1.01*tnu.E2nu || rfield.anu[nu-1]<tnu.E2nu );
00100
00101
00102
00103 sum += tnu.As2nu[nu];
00104
00105
00106
00107 if( lgDoInduced )
00108 {
00109
00110 double rate_up = 0.5 * tnu.As2nu[nu] *
00111 rfield.SummedOcc[nu] * rfield.SummedOcc[tnu.ipSym2nu[nu]-1];
00112 tnu.induc_up += rate_up;
00113 tnu.induc_dn += rate_up + tnu.As2nu[nu] *
00114 (rfield.SummedOcc[nu] + rfield.SummedOcc[tnu.ipSym2nu[nu]-1]);
00115 }
00116 }
00117
00118
00119
00120 ASSERT( fabs( 1.f - (realnum)sum/tnu.AulTotal ) < 0.01f );
00121
00122 return;
00123 }
00124
00125 void CalcTwoPhotonEmission( two_photon& tnu, bool lgDoInduced )
00126 {
00127 DEBUG_ENTRY( "CalcTwoPhotonEmission()" );
00128
00129
00130 ASSERT( tnu.ipTwoPhoE>0 );
00131
00132
00133
00134 for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
00135 {
00136
00137
00138 tnu.local_emis[nu] = 2.f * (realnum)(*tnu.Pop) * tnu.As2nu[nu];
00139 }
00140
00141
00142
00143 if( lgDoInduced )
00144 {
00145 for( long nu=0; nu < tnu.ipTwoPhoE; nu++ )
00146 {
00147
00148 tnu.local_emis[nu] *= (1.f + rfield.SummedOcc[nu]) *
00149 (1.f+rfield.SummedOcc[tnu.ipSym2nu[nu]-1]);
00150 }
00151 }
00152
00153 return;
00154 }
00155
00156
00157 void PrtTwoPhotonEmissCoef( const two_photon& tnu, const double& densityProduct )
00158 {
00159 DEBUG_ENTRY( "PrtTwoPhotonEmissCoef()" );
00160
00161 fprintf( ioQQQ, "\ny\tGammaNot(2q)\n");
00162
00163 for( long yTimes20=1; yTimes20<=10; yTimes20++ )
00164 {
00165 double y = yTimes20/20.;
00166
00167 fprintf( ioQQQ, "%.3e\t", (double)y );
00168
00169 long i = ipoint(y*tnu.E2nu);
00170 fprintf( ioQQQ, "%.3e\n",
00171 8./3.*HPLANCK*(*tnu.Pop)/densityProduct*y*tnu.As2nu[i]*tnu.E2nu/rfield.widflx[i] );
00172 }
00173
00174 return;
00175 }
00176