Changeset 2112 for branches/newmole

Timestamp:

06/17/08 17:10:09 (7 months ago)

Author:

rjrw

Message:

Merged from trunk r2088:2111

Location:

branches/newmole/source

Files:

• atmdat_adfa.cpp (modified) (1 diff)
• cddefines.h (modified) (2 diffs)
• container_classes.h (modified) (16 diffs)
• grains.cpp (modified) (14 diffs)
• grains_mie.cpp (modified) (3 diffs)
• grains_qheat.cpp (modified) (10 diffs)
• grainvar.h (modified) (1 diff)
• ion_solver.cpp (modified) (2 diffs)
• iso.h (modified) (1 diff)
• iso_ionize_recombine.cpp (modified) (1 diff)
• iso_level.cpp (modified) (3 diffs)
• Makefile (modified) (6 diffs)
• mole_species.cpp (modified) (1 diff)
• parse_grid.cpp (modified) (1 diff)
• parse_punch.cpp (modified) (1 diff)
• punch_do.cpp (modified) (2 diffs)

branches/newmole/source/atmdat_adfa.cpp

@@ -236 +236 @@
         ASSERT( !lgErr );
 
+        /*refer HI      cs      Anderson, H., Ballance, C.P., Badnell, N.R., 
+         *refercon      & Summers, H.P  2000, J Phys B, 33, 1255 */
         const static char chFile8[] = "h_coll_str.dat";
 

branches/newmole/source/cddefines.h

@@ -69 +69 @@
 typedef float realnum;
 #endif
+
+//Compile-time assertion after Alexandrescu
+template<bool> struct StaticAssertFailed;
+template<> struct StaticAssertFailed<true> {};
+#define STATIC_ASSERT(x) (StaticAssertFailed< (x) == true >())
 
 typedef float sys_float;
@@ -1338 +1343 @@
 
 /* Explicit instantiations for debugging purposes */
-INSTANTIATE_MULTI_ARR( quantumState, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL );
-INSTANTIATE_MULTI_ARR( transition, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL );
+INSTANTIATE_MULTI_ARR( quantumState, lgBOUNDSCHECKVAL );
+INSTANTIATE_MULTI_ARR( transition, lgBOUNDSCHECKVAL );
 
 

branches/newmole/source/container_classes.h

@@ -32 +32 @@
 //! basic_pntr - base class for generalization of normal pointers that enables bounds checking
 //! it comes with the full set of operators that you would expect for a random access pointer
-template<class T, int d, mem_layout ALLOC, bool lgBC>
+template<class T, bool lgBC>
 class basic_pntr
 {
@@ -126 +126 @@
 
 //! pntr - interface class to replace normal pointers
-template<class T, int d, mem_layout ALLOC, bool lgBC>
-class pntr : public basic_pntr<T,d,ALLOC,lgBC>
+template<class T, bool lgBC>
+class pntr : public basic_pntr<T,lgBC>
 {
 public:
         // constructors are not inherited, so define them again
-        pntr( T* p0 ) : basic_pntr<T,d,ALLOC,lgBC>( p0 ) {}
-        pntr( T* p0, T* p1, T* p2 ) : basic_pntr<T,d,ALLOC,lgBC>( p0, p1, p2 ) {}
+        pntr( T* p0 ) : basic_pntr<T,lgBC>( p0 ) {}
+        pntr( T* p0, T* p1, T* p2 ) : basic_pntr<T,lgBC>( p0, p1, p2 ) {}
         pntr() {}
         // the increment / decrement operators need to be recast...
         // otherwise expressions like p = ++q would be illegal for iterators...
-        pntr& operator++ () { return static_cast<pntr&>(basic_pntr<T,d,ALLOC,lgBC>::operator++()); }
+        pntr& operator++ () { return static_cast<pntr&>(basic_pntr<T,lgBC>::operator++()); }
         const pntr operator++ (int) { pntr t = *this; ++(*this); return t; }
-        pntr& operator-- () { return static_cast<pntr&>(basic_pntr<T,d,ALLOC,lgBC>::operator--()); }
+        pntr& operator-- () { return static_cast<pntr&>(basic_pntr<T,lgBC>::operator--()); }
         const pntr operator-- (int) { pntr t = *this; --(*this); return t; }
         // define p+n, p-n, p-q
@@ -147 +147 @@
 
 // this defines n+p
-template<class T, int d, mem_layout ALLOC, bool lgBC>
-inline const pntr<T,d,ALLOC,lgBC> operator+ ( const ptrdiff_t n, const pntr<T,d,ALLOC,lgBC>& t )
-{
-        pntr<T,d,ALLOC,lgBC> s = t;
+template<class T, bool lgBC>
+inline const pntr<T,lgBC> operator+ ( const ptrdiff_t n, const pntr<T,lgBC>& t )
+{
+        pntr<T,lgBC> s = t;
         s += n;
         return s;
@@ -156 +156 @@
 
 //! const_pntr - same as pntr, except that it replaces const pointers rather than normal pointers
-template<class T, int d, mem_layout ALLOC, bool lgBC>
-class const_pntr : public basic_pntr<T,d,ALLOC,lgBC>
+template<class T, bool lgBC>
+class const_pntr : public basic_pntr<T,lgBC>
 {
 public:
         // constructors are not inherited, so define them again
-        const_pntr( T* p0 ) : basic_pntr<T,d,ALLOC,lgBC>( p0 ) {}
-        const_pntr( T* p0, T* p1, T* p2 ) : basic_pntr<T,d,ALLOC,lgBC>( p0, p1, p2 ) {}
+        const_pntr( T* p0 ) : basic_pntr<T,lgBC>( p0 ) {}
+        const_pntr( T* p0, T* p1, T* p2 ) : basic_pntr<T,lgBC>( p0, p1, p2 ) {}
         const_pntr() {}
         // make sure we can assign a pntr to a const_pntr by creating an implicit conversion to const_pntr
-        const_pntr( const pntr<T,d,ALLOC,lgBC>& t ) : basic_pntr<T,d,ALLOC,lgBC>( t ) {}
+        const_pntr( const pntr<T,lgBC>& t ) : basic_pntr<T,lgBC>( t ) {}
         // the increment / decrement operators need to be recast...
         // otherwise expressions like *p++ = 1. would be legal for const_iterators...
-        const_pntr& operator++ () { return static_cast<const_pntr&>(basic_pntr<T,d,ALLOC,lgBC>::operator++()); }
+        const_pntr& operator++ () { return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator++()); }
         const const_pntr operator++ (int) { const_pntr t = *this; ++(*this); return t; }
-        const_pntr& operator-- () { return static_cast<const_pntr&>(basic_pntr<T,d,ALLOC,lgBC>::operator--()); }
+        const_pntr& operator-- () { return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator--()); }
         const const_pntr operator-- (int) { const_pntr t = *this; --(*this); return t; }
         const_pntr& operator+= ( const ptrdiff_t n )
         {
-                return static_cast<const_pntr&>(basic_pntr<T,d,ALLOC,lgBC>::operator+=(n));
+                return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator+=(n));
         }
         const_pntr& operator-= ( const ptrdiff_t n )
         {
-                return static_cast<const_pntr&>(basic_pntr<T,d,ALLOC,lgBC>::operator-=(n));
+                return static_cast<const_pntr&>(basic_pntr<T,lgBC>::operator-=(n));
         }
         // the dereference operators need to be recast...
-        const T& operator* () const { return static_cast<const T&>(basic_pntr<T,d,ALLOC,lgBC>::operator*()); }
-        const T* operator-> () const { return static_cast<const T*>(basic_pntr<T,d,ALLOC,lgBC>::operator->()); }
+        const T& operator* () const { return static_cast<const T&>(basic_pntr<T,lgBC>::operator*()); }
+        const T* operator-> () const { return static_cast<const T*>(basic_pntr<T,lgBC>::operator->()); }
         const T& operator[] ( const ptrdiff_t n ) const
         {
-                return static_cast<const T&>(basic_pntr<T,d,ALLOC,lgBC>::operator[](n));
+                return static_cast<const T&>(basic_pntr<T,lgBC>::operator[](n));
         }
         // define p+n, p-n, p-q
@@ -194 +194 @@
 
 // this defines n+p
-template<class T, int d, mem_layout ALLOC, bool lgBC>
-inline const const_pntr<T,d,ALLOC,lgBC> operator+ ( const ptrdiff_t n, const const_pntr<T,d,ALLOC,lgBC>& t )
-{
-        const_pntr<T,d,ALLOC,lgBC> s = t;
+template<class T, bool lgBC>
+inline const const_pntr<T,lgBC> operator+ ( const ptrdiff_t n, const const_pntr<T,lgBC>& t )
+{
+        const_pntr<T,lgBC> s = t;
         s += n;
         return s;
@@ -275 +275 @@
 //! multi_geom - this class maintains all the geometry information for multi_arr
 //! keeping it separate makes it easy to clone the information from one multi_arr to another
-template<int d,mem_layout ALLOC=MEM_LAYOUT_VAL,bool lgBC=lgBOUNDSCHECKVAL>
+template<int d,mem_layout ALLOC=MEM_LAYOUT_VAL>
 class multi_geom
 {
@@ -382 +382 @@
 #       pragma warning( disable : 4127 )
 #endif
+
+                STATIC_ASSERT ( ALLOC == C_TYPE || ALLOC == ARPA_TYPE );
+
                 if( ALLOC == ARPA_TYPE )
                 {
@@ -403 +406 @@
                 {
                         TotalInsanity();
-                }               
+                }
         }
 
@@ -902 +905 @@
 {
         // ancillary data describing the memory layout of the multi_arr
-        multi_geom<d,ALLOC,lgBC> p_g;
+        multi_geom<d,ALLOC> p_g;
         T** p_psl[d-1];     // pointer arrays for ARPA structure
         valarray<T> p_dsl;  // this contains the actual data
@@ -921 +924 @@
         typedef size_t       size_type;
         typedef ptrdiff_t    difference_type;
-        typedef pntr<T,d,ALLOC,lgBC>       iterator;
-        typedef const_pntr<T,d,ALLOC,lgBC> const_iterator;
+        typedef pntr<T,lgBC>       iterator;
+        typedef const_pntr<T,lgBC> const_iterator;
 
 private:
@@ -951 +954 @@
                 p_clear1();
         }
-        multi_arr(const multi_geom<d,ALLOC,lgBC>& g)
+        multi_arr(const multi_geom<d,ALLOC>& g)
         {
                 p_clear1();
@@ -1026 +1029 @@
 #       pragma warning( disable : 4127 )
 #endif
+                STATIC_ASSERT ( ALLOC == C_TYPE || ALLOC == ARPA_TYPE );
+
                 if( ALLOC == ARPA_TYPE )
                 {
@@ -1064 +1069 @@
         }
         // clone the geometry from another multi_arr
-        void alloc(const multi_geom<d,ALLOC,lgBC>& g)
+        void alloc(const multi_geom<d,ALLOC>& g)
         {
                 if( &g != &p_g )
@@ -1647 +1652 @@
         }
 
-        const multi_geom<d,ALLOC,lgBC>& clone() const
+        const multi_geom<d,ALLOC>& clone() const
         {
                 return p_g;
@@ -1719 +1724 @@
 
 // on Mac systems these instantiations need to be extern in order to avoid duplicate symbols
-#define INSTANTIATE_MULTI_ARR( TYPE, LAYOUT, BC ) \
-INST_EXTERN template class pntr<TYPE,2,LAYOUT,BC>; \
-INST_EXTERN template class pntr<TYPE,3,LAYOUT,BC>; \
-INST_EXTERN template class pntr<TYPE,4,LAYOUT,BC>; \
-INST_EXTERN template class pntr<TYPE,5,LAYOUT,BC>; \
-INST_EXTERN template class pntr<TYPE,6,LAYOUT,BC>; \
-INST_EXTERN template class const_pntr<TYPE,2,LAYOUT,BC>; \
-INST_EXTERN template class const_pntr<TYPE,3,LAYOUT,BC>; \
-INST_EXTERN template class const_pntr<TYPE,4,LAYOUT,BC>; \
-INST_EXTERN template class const_pntr<TYPE,5,LAYOUT,BC>; \
-INST_EXTERN template class const_pntr<TYPE,6,LAYOUT,BC>
-
-INSTANTIATE_MULTI_ARR( bool, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL );
-INSTANTIATE_MULTI_ARR( long, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL );
-INSTANTIATE_MULTI_ARR( realnum, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL );
-INSTANTIATE_MULTI_ARR( double, MEM_LAYOUT_VAL, lgBOUNDSCHECKVAL );
-INSTANTIATE_MULTI_ARR( double, C_TYPE, lgBOUNDSCHECKVAL );
+#define INSTANTIATE_MULTI_ARR( TYPE, BC ) \
+INST_EXTERN template class pntr<TYPE,BC>; \
+INST_EXTERN template class const_pntr<TYPE,BC>;
+
+INSTANTIATE_MULTI_ARR( bool, lgBOUNDSCHECKVAL );
+INSTANTIATE_MULTI_ARR( long, lgBOUNDSCHECKVAL );
+INSTANTIATE_MULTI_ARR( realnum,lgBOUNDSCHECKVAL );
+INSTANTIATE_MULTI_ARR( double, lgBOUNDSCHECKVAL );
 
 
@@ -1758 +1754 @@
         typedef long         size_type;
         typedef ptrdiff_t    difference_type;
-        typedef pntr<T,1,FLX_TYPE,lgBC>       iterator;
-        typedef const_pntr<T,1,FLX_TYPE,lgBC> const_iterator;
+        typedef pntr<T,lgBC>       iterator;
+        typedef const_pntr<T,lgBC> const_iterator;
 
 private:

branches/newmole/source/grains.cpp

@@ -33 +33 @@
 inline double ASINH(double x)
 {
-        if( abs(x) <= 2.3e-4 )
-                return x - pow3(x)/6.;
-        else if( abs(x) <= 8.e-3 )
+        if( abs(x) <= 8.e-3 )
                 return ((0.075*pow2(x) - 1./6.)*pow2(x) + 1.)*x;
         else if( abs(x) <= 1./sqrt(DBL_EPSILON) )
@@ -47 +45 @@
         {
                 if( x < 0. )
-                        return -log(-2.*x);
+                        return -(log(-x)+LN_TWO);
                 else
-                        return log(2.*x);
+                        return log(x)+LN_TWO;
         }
 }
@@ -145 +143 @@
                 return 1.e-7;
         else
-                return 3.e-6*pow((double)gv.bin[nd]->dustp[0],-0.85)*pow(e*EVRYD*1.e-3,1.5);
+                return 3.e-6*gv.bin[nd]->eec*sqrt(pow3(e*EVRYD*1.e-3));
 }
         
@@ -192 +190 @@
 STATIC void UpdatePot2(long,long);
 /* Helper function to calculate primary and secondary yields and the average electron energy at infinity */
-STATIC void Yfunc(long,long,double,double,double,double,double,/*@out@*/double*,/*@out@*/double*,
+inline void Yfunc(long,long,double,double,double,double,double,/*@out@*/double*,/*@out@*/double*,
                   /*@out@*/double*,/*@out@*/double*);
 /* This calculates the y0 function for band electrons (Sect. 4.1.3/4.1.4 of WDB06) */
@@ -203 +201 @@
 STATIC double y1psa(long,long,double);
 /* This calculates the y2 function for primary and Auger electrons (Eq. 8 of WDB06) */
-STATIC double y2pa(double,double,long,double*);
+inline double y2pa(double,double,long,double*);
 /* This calculates the y2 function for secondary electrons (Eq. 20-21 of WDB06) */
-STATIC double y2s(double,double,long,double*);
+inline double y2s(double,double,long,double*);
 /* find highest ionization stage with non-zero population */
 STATIC long HighestIonStage(void);
@@ -3216 +3214 @@
         }
 
+        double GrainPot = chrg2pot(Zg,nd);
+
         if( nfill > ipLo )
         {
@@ -3241 +3241 @@
                         /* calculate yield for band structure */
                         Ehi = Ehi_band = anu[i] - ThresInfVal - Emin;
-                        Wcorr = ThresInfVal + Emin - chrg2pot(Zg,nd);
+                        Wcorr = ThresInfVal + Emin - GrainPot;
                         Eel = anu[i] - DustWorkFcn;
                         yzero = y0b( nd, nz, i );
@@ -3252 +3252 @@
 
                         /* add in yields for inner shell ionization */
-                        for( ns=1; ns < gv.bin[nd]->nShells; ns++ )
+                        unsigned long nsmax = gv.bin[nd]->nShells;
+                        for( ns=1; ns < nsmax; ns++ )
                         {
                                 sptr = gv.bin[nd]->sd[ns];
@@ -3268 +3269 @@
 
                                 /* add in Auger yields */
-                                for( n=0; n < sptr->nData; n++ )
+                                long nmax = sptr->nData;
+                                for( n=0; n < nmax; n++ )
                                 {
                                         double max = sptr->AvNr[n]*sptr->p[i];
-                                        Ehi = sptr->Ener[n] - chrg2pot(Zg,nd);
+                                        Ehi = sptr->Ener[n] - GrainPot;
                                         Eel = sptr->Ener[n];
                                         Yfunc(nd,nz,sptr->y01A[n][i],max,Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
@@ -3432 +3434 @@
 
 /* Helper function to calculate primary and secondary yields and the average electron energy at infinity */
-STATIC void Yfunc(long nd,
+inline void Yfunc(long nd,
                   long nz,
                   double y01,
@@ -3472 +3474 @@
                 /* this is Eq. 18 of WDB06 */
                 /* Eel may be negative near threshold -> set yield to zero */
-                f3 = max(Eel,0.)/(eps*elec_esc_length(Eel,nd)*(1./gv.bin[nd]->AvRadius + 1.e7));
+                f3 = max(Eel,0.)/(eps*elec_esc_length(Eel,nd)*gv.bin[nd]->eyc);
                 *Ys = y2sec*f3*min(y01,maxval);
         }
@@ -3614 +3616 @@
 
 /* This calculates the y2 function for primary and Auger electrons (Eq. 8 of WDB06) */
-STATIC double y2pa(double Elo,
+inline double y2pa(double Elo,
                    double Ehi,
                    long Zg,
@@ -3658 +3660 @@
 
 /* This calculates the y2 function for secondary electrons (Eqs. 20-21 of WDB06) */
-STATIC double y2s(double Elo,
+inline double y2s(double Elo,
                   double Ehi,
                   long Zg,
@@ -3713 +3715 @@
                                 {
                                         N0 = alpha*(1./sqR0 - 1./sqRh);
-                                        E0 = alpha*ETILDE*(ASINH(x) + ASINH(yl) - yh/sqRh);
+                                        E0 = alpha*ETILDE*(ASINH(x*sqR0 + yl*sqRh) - yh/sqRh);
                                 }
                         }

branches/newmole/source/grains_mie.cpp

@@ -1088 +1088 @@
         mie_next_data(chFile,io2,chLine,&dl);
         mie_read_realnum(chFile,chLine,&gv.bin[nd]->dustp[0],true,dl);
+        /* constant needed in the evaluation of the electron escape length */
+        gv.bin[nd]->eec = pow((double)gv.bin[nd]->dustp[0],-0.85);
 
         /* molecular weight of grain molecule (amu) */
@@ -1108 +1110 @@
         mie_next_data(chFile,io2,chLine,&dl);
         mie_read_realnum(chFile,chLine,&gv.bin[nd]->AvRadius,true,dl);
+        gv.bin[nd]->eyc = 1./gv.bin[nd]->AvRadius + 1.e7;
 
         /* average grain area <4pi*a^2> for entire size distr (cm^2) */
@@ -1279 +1282 @@
                                         gv.bin[nd2]->dustp[j] = gv.bin[nd]->dustp[j];
                                 }
+                                gv.bin[nd2]->eec = gv.bin[nd]->eec;
+                                gv.bin[nd2]->eyc = gv.bin[nd]->eyc;
                                 gv.bin[nd2]->DustWorkFcn = gv.bin[nd]->DustWorkFcn;
                                 gv.bin[nd2]->BandGap = gv.bin[nd]->BandGap;

branches/newmole/source/grains_qheat.cpp

@@ -622 +622 @@
         if( trace.lgTrace && trace.lgDustBug )
         {
-                fprintf( ioQQQ, "   grain heating: %.4e, integral %.4e, total rate %.4e\n",
-                         gv.bin[nd]->GrainHeat,integral,Phi[0]);
+                double Rate2 = 0.;
+                for( int nz=0; nz < gv.bin[nd]->nChrg; nz++ )
+                        Rate2 += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->HeatingRate2;
+
+                fprintf( ioQQQ, "   grain heating: %.4e, integral %.4e, total rate %.4e lgNegRate %c\n",
+                         gv.bin[nd]->GrainHeat,integral,Phi[0],TorF(lgNegRate));
                 fprintf( ioQQQ, "   av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
                          gv.bin[nd]->tedust,Umax);
                 fprintf( ioQQQ, "   fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
                          NumEvents,fwhm,FwhmRatio );
-                fprintf( ioQQQ, "   lgQHTooWide %d\n", gv.bin[nd]->lgQHTooWide );
+                fprintf( ioQQQ, "   HeatingRate1 %.4e HeatingRate2 %.4e lgQHTooWide %c\n",
+                         gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3, Rate2*gv.bin[nd]->cnv_H_pCM3,
+                         TorF(gv.bin[nd]->lgQHTooWide) );
         }
 
@@ -911 +917 @@
          * energies, the reshuffling for higher energies is done in the next loop
          * phiTilde has units events/H/s/cell at default depletion */
+
+        double NegHeatingRate = 0.;
 
         for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
@@ -983 +991 @@
                 *check += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
 
+                sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
+
                 if( DEBUG_LOC )
                 {
@@ -990 +1000 @@
                                 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
                         }
-                        sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
                         dprintf( ioQQQ, " integral test 1: integral %.6e %.6e\n", integral, sum );
                 }
@@ -1007 +1016 @@
 
                 /* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, pah_crash.in, PvH */
-                if( fabs(gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3) > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat ) 
+                if( gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat ) 
                 {
                         double *RateArr,Sum,ESum,DSum,Delta,E_av2,Corr;
@@ -1076 +1085 @@
                         }
 
+                        sum += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
+
                         if( DEBUG_LOC )
                         {
@@ -1083 +1094 @@
                                         integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
                                 }
-                                sum += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
                                 dprintf( ioQQQ, " integral test 2: integral %.6e %.6e\n", integral, sum );
                         }
 
                         free( RateArr );
+                }
+                else
+                {
+                        NegHeatingRate += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
                 }
         }
@@ -1109 +1123 @@
 
         /* >>chng 03 nov 06, check for extremely low HeatingRate and save CPU time, PvH */
-        if( fabs(gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3) > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
+        if( gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
         {
                 /* limits for Taylor expansion of (1+x)*exp(-x) */
@@ -1153 +1167 @@
                         ylo = yhi;
                 }
-        }
-
-        if( DEBUG_LOC )
-        {
-                double integral = 0.;
+
+                sum += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
+
+                if( DEBUG_LOC )
+                {
+                        double integral = 0.;
+                        for( i=0; i < gv.bin[nd]->qnflux; i++ )
+                        {
+                                integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
+                        }
+                        dprintf( ioQQQ, " integral test 3: integral %.6e %.6e\n", integral, sum );
+                }
+        }
+        else
+        {
+                NegHeatingRate += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
+        }
+
+        /* here we account for the negative heating rates, we simply do that by scaling the entire
+         * phiTilde array down by a constant factor such that the total amount of energy is conserved
+         * This treatment assures that phiTilde never goes negative, which avoids problems further on */
+        if( NegHeatingRate < 0. )
+        {
+                double scale_fac = (sum+NegHeatingRate)/sum;
                 for( i=0; i < gv.bin[nd]->qnflux; i++ )
-                {
-                        integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
-                }
-                sum += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
-                dprintf( ioQQQ, " integral test 3: integral %.6e %.6e\n", integral, sum );
-        }
+                        phiTilde[i] *= scale_fac;
+
+                sum += NegHeatingRate;
+
+                if( DEBUG_LOC )
+                {
+                        double integral = 0.;
+                        for( i=0; i < gv.bin[nd]->qnflux; i++ )
+                        {
+                                integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
+                        }