Changeset 2151

Timestamp:

07/03/08 10:15:35 (8 weeks ago)

Author:

gary

Message:

bugfix - factor of 1/speedlight /xMassDensity was missing from force multiplier - this is only used in some print and punch statements so did not affect the calculation.

Files:

• trunk/source/radius_increment.cpp (modified) (11 diffs)

trunk/source/radius_increment.cpp

@@ -186 +186 @@
                 {
                         fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], 
-                          rfield.flux[i] + rfield.outlin[i] + rfield.ConInterOut[i], 
+                          rfield.flux[0][i] + rfield.outlin[0][i] + rfield.ConInterOut[i], 
                           rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]);
                 }
@@ -202 +202 @@
         for( i=0; i < rfield.nflux; i++ )
         {
-                if( rfield.flux[i] < 0. )
+                if( rfield.flux[0][i] < 0. )
                 {
                         fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
                         fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 
-                          rfield.flux[i], rfield.anu[i], i );
+                          rfield.flux[0][i], rfield.anu[i], i );
                         lgFlxNeg = true;
                 }
@@ -230 +230 @@
                         lgFlxNeg = true;
                 }
-                if( rfield.outlin[i] < 0. )
+                if( rfield.outlin[0][i] < 0. )
                 {
                         fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
                         fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 
-                          rfield.outlin[i], rfield.anu[i], i, rfield.chLineLabel[i]  );
+                          rfield.outlin[0][i], rfield.anu[i], i, rfield.chLineLabel[i]  );
                         lgFlxNeg = true;
                 }
@@ -346 +346 @@
                  * drNext was set by NEXTDR and will be next dr */
                 /* compute both total and Thomson scat rad acceleration */
-                rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ 
+                rforce += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+ 
                         rfield.ConInterOut[i])* rfield.anu[i]*(opac.opacity_abs[i] + 
                   opac.opacity_sct[i]);
 
-                relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ 
+                relec += ((rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+ 
                         rfield.ConInterOut[i])* escatc)*rfield.anu[i];
 
@@ -383 +383 @@
                 rfield.flux_beam_time[i] *= (realnum)opfac;
                 rfield.flux_isotropic[i] *= (realnum)opfac;
-                rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
+                rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
                         rfield.flux_isotropic[i];
 
                 /* update SummedCon here since flux changed */
-                rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
+                rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
 
                 /* outward lines and continua */
                 rfield.ConInterOut[i] *= (realnum)opfac;
-                rfield.outlin[i] *= (realnum)opfac;
+                rfield.outlin[0][i] *= (realnum)opfac;
                 rfield.outlin_noplot[i] *= (realnum)opfac;
 
-                rfield.ConEmitOut[i] *= (realnum)opfac;
-                rfield.ConEmitOut[i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i]/**/;
+                rfield.ConEmitOut[0][i] *= (realnum)opfac;
+                rfield.ConEmitOut[0][i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i]/**/;
 
                 /* set occupation numbers, first attenuated incident continuum */
-                rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i];
+                rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
 
                 /* >>chng 00 oct 03, add diffuse continua */
@@ -405 +405 @@
 
                 /* >>chng 05 feb 16, occupation number of outward diffuse continuum */
-                rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[i]*rfield.convoc[i];
+                rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[0][i]*rfield.convoc[i];
         }
 
@@ -416 +416 @@
         /* >>chng 02 dec 05, add test for small float, protect against models 
          * where we have gone below smallfloat, and so float is ragged */
-        if( rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT &&
-                (rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
-                SDIV(rfield.flux_total_incident[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) > SMALLFLOAT && 
+        if( rfield.flux[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT &&
+                (rfield.flux[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
+                SDIV(rfield.flux_total_incident[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) > SMALLFLOAT && 
                 !opac.lgScatON &&
                 radius.depth/radius.Radius < 1e-4 )
@@ -427 +427 @@
                  * and here convert to double.
                  * error found by Peter van Hoof */
-                double tau_effec = -log((double)rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
+                double tau_effec = -log((double)rfield.flux[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
                         (double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
-                        (double)rfield.flux_total_incident[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]);
+                        (double)rfield.flux_total_incident[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]);
 
                 /* this is computed absorption optical depth to illuminated face */
@@ -472 +472 @@
                         rfield.flux_beam_time[i] *= (realnum)opfac;
                         rfield.flux_isotropic[i] *= (realnum)opfac;
-                        rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
+                        rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
                                 rfield.flux_isotropic[i];
 
                         rfield.ConInterOut[i] *= (realnum)opfac;
-                        rfield.ConEmitOut[i] *= (realnum)opfac;
-                        rfield.outlin[i] *= (realnum)opfac;
+                        rfield.ConEmitOut[0][i] *= (realnum)opfac;
+                        rfield.outlin[0][i] *= (realnum)opfac;
                         rfield.outlin_noplot[i] *= (realnum)opfac;
                 }
@@ -485 +485 @@
         cmshft();
 
-        relec *= EN1RYD;
-
-#       if 0
-        /* this is all evaluated in pressure_total */
-        /* radiative acceleration; xMassDensity is gm per cc, evaluated when PTOT called */
-        t = 1./SPEEDLIGHT/dense.xMassDensity;
-        /* >>chng 03 feb 01, lgLineRadPresOn was missing */
-        wind.AccelLine = (realnum)(RT_line_driving()*t)*pressure.lgLineRadPresOn;
-        /* >>chng 01 aug 03, add lgConPress to kill AccelCont here */
-        wind.AccelCont = (realnum)(rforce*EN1RYD*t)*pressure.lgContRadPresOn;
-        /* this is numerically unstable */
-        wind.AccelPres = 0.;
-        /* total acceleration */
-        wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres;
-#       endif
         /* remember the largest value */
         wind.AccelMax = (realnum)MAX2(wind.AccelMax,wind.AccelTot);
 
-        /* force multiplier; RELEC*T can be zero for very low densities */
-        if( relec*t > SMALLFLOAT )
-        {
-                wind.fmul = (realnum)((wind.AccelLine + wind.AccelCont)/(relec*t));
-        }
-        else
-        {
-                wind.fmul = 0.;
-        }
+        /* force multiplier; relec can be zero for very low densities - so use double
+         * form of safe_div - fmul is of order unity - wind.AccelLine and wind.AccelCont 
+         * are both floats to will underflow long before relec will - fmul is only use
+         * in output, not in any physics */
+        relec *= (EN1RYD/SPEEDLIGHT/dense.xMassDensity);
+        wind.fmul = (realnum)safe_div( (double)(wind.AccelLine + wind.AccelCont) , 
+                relec );
 
         /* keep track of average acceleration */
@@ -558 +541 @@
 
                                 /* ConEmitReflec - reflected diffuse continuum */
-                                rfield.ConEmitReflec[i] += (realnum)(Reflec_Diffuse_Cont);
+                                rfield.ConEmitReflec[0][i] += (realnum)(Reflec_Diffuse_Cont);
 
                                 /* the reflected part of the incident continuum */
-                                rfield.ConRefIncid[i] += (realnum)(rfield.flux[i]*opac.opacity_sct[i]*
+                                rfield.ConRefIncid[0][i] += (realnum)(rfield.flux[0][i]*opac.opacity_sct[i]*
                                   radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
 
                                 /* reflected line emission */
-                                rfield.reflin[i] += (realnum)(rfield.outlin[i]*opac.opacity_sct[i]*
+                                rfield.reflin[0][i] += (realnum)(rfield.outlin[0][i]*opac.opacity_sct[i]*
                                   radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
                         }