cool_save.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*CoolSave save coolants */
#include "cddefines.h"
#include "thermal.h"
#include "dynamics.h"
#include "radius.h"
#include "conv.h"
#include "phycon.h"
#include "save.h"
#include "grainvar.h"
#include "hmi.h"
#include "coolheavy.h"
 
 
/* this is limit to number of coolants to print out */
static const int IPRINT = 100;
 
/*CoolSave save coolants */
void CoolSave(FILE * io, char chJob[])
{
        long int i, 
          ip, 
          is;
 
        int nFail;
 
        double cset, 
                cool_total,
                heat_total;
 
        realnum
                *csav,
                *sgnsav;
        long int *index;
 
        DEBUG_ENTRY( "CoolSave()" );
 
        /* cannot do one-time init since thermal.ncltot can change */
        index = (long int *)CALLOC((size_t)thermal.ncltot,sizeof(long int));
        csav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum));
        sgnsav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum));
 
        cool_total = thermal.ctot;
        heat_total = thermal.htot;
 
        /* >>chng 06 mar 17, comment out following block and replace with this 
         * removing dynamics heating & cooling and report only physical
         * heating and cooling 
         * NB the heating and cooling as punched no longer need be
         * equal for a converged model */
        cool_total -= dynamics.Cool();
        heat_total -= dynamics.Heat();
#       if 0
        if(dynamics.Cool > dynamics.Heat()) 
        {
                cool_total -= dynamics.Heat();
                heat_total -= dynamics.Heat();
        } 
        else
        {
                cool_total -= dynamics.Cool;
                heat_total -= dynamics.Cool;
        }
#       endif
 
        /* cset will be weakest cooling to consider
         * WeakHeatCool set with 'set weakheatcool' command
         * default is 0.05 */
        cset = cool_total*save.WeakHeatCool;
 
        /* first find all strong lines, both + and - sign */
        ip = thermal.ncltot;
 
        for( i=0; i < ip; i++ )
        {
                csav[i] = (realnum)(MAX2(thermal.cooling[i],thermal.heatnt[i])/
                        SDIV(cool_total));
 
                /* save sign to remember if heating or cooling line */
                if( thermal.heatnt[i] == 0. )
                {
                        sgnsav[i] = 1.;
                }
                else
                {
                        sgnsav[i] = -1.;
                }
        }
 
        /* order strongest to weakest */
        /* now sort by decreasing importance */
        /*spsort netlib routine to sort array returning sorted indices */
        spsort(
                  /* input array to be sorted */
                  csav, 
                  /* number of values in x */
                  ip, 
                  /* permutation output array */
                  index, 
                  /* flag saying what to do - 1 sorts into increasing order, not changing
                   * the original routine */
                  -1, 
                  /* error condition, should be 0 */
                  &nFail);
 
        /* warn if tcovergence failure occurred */
        if( !conv.lgConvTemp )
        {
                fprintf( io, "#>>>>  Temperature not converged.\n" );
        }
        else if( !conv.lgConvEden )
        {
                fprintf( io, "#>>>>  Electron density not converged.\n" );
        }
        else if( !conv.lgConvIoniz() )
        {
                fprintf( io, "#>>>>  Ionization not converged.\n" );
        }
        else if( !conv.lgConvPres )
        {
                fprintf( io, "#>>>>  Pressure not converged.\n" );
        }
 
        if( strcmp(chJob,"EACH") == 0 )
        {
                /* begin the print out with zone number, total heating and cooling */
                fprintf( io, "%.5e\t%.4e\t%.4e", 
                                 radius.depth_mid_zone, 
                                 phycon.te, 
                                 cool_total );
                double debug_ctot = 0.;
                
                for( int i=0 ; i <= LIMELM ; i++)
                {
                        fprintf( io, "\t%.4e", thermal.elementcool[i] );
                        debug_ctot += thermal.elementcool[i];
                }
                
                /*fprintf( io, "\t%.4e", thermal.dust );
                fprintf( io, "\t%.4e", thermal.H2cX );
                fprintf( io, "\t%.4e", thermal.CT_C );
                fprintf( io, "\t%.4e", thermal.H_fb );
                fprintf( io, "\t%.4e", thermal.H2ln );
                fprintf( io, "\t%.4e", thermal.HDro );
                fprintf( io, "\t%.4e", thermal.H2p );
                fprintf( io, "\t%.4e", thermal.FF_c );
                fprintf( io, "\t%.4e", thermal.hvFB );
                fprintf( io, "\t%.4e", thermal.eeff );
                fprintf( io, "\t%.4e", thermal.adve );
                fprintf( io, "\t%.4e", thermal.Comp );
                fprintf( io, "\t%.4e", thermal.Extr );
                fprintf( io, "\t%.4e", thermal.Expn );
                fprintf( io, "\t%.4e", thermal.Cycl );
                fprintf( io, "\t%.4e", thermal.Hvin );
                fprintf( io, " \n" );*/
 
                fprintf( io, "\t%.4e", MAX2(0.,gv.GasCoolColl) );
                debug_ctot += MAX2(0.,gv.GasCoolColl);
                
                fprintf( io, "\t%.4e", MAX2(0.,-hmi.HeatH2Dexc_used) );
                debug_ctot += MAX2(0.,-hmi.HeatH2Dexc_used);
                
                fprintf( io, "\t%.4e", thermal.char_tran_cool );
                debug_ctot += thermal.char_tran_cool;
                
                fprintf( io, "\t%.4e", hmi.hmicol );
                debug_ctot += hmi.hmicol;
                
                fprintf( io, "\t%.4e", CoolHeavy.h2line );
                debug_ctot += CoolHeavy.h2line;
                
                fprintf( io, "\t%.4e", CoolHeavy.HD );
                debug_ctot += CoolHeavy.HD;
                
                fprintf( io, "\t%.4e", CoolHeavy.H2PlsCool );
                debug_ctot += CoolHeavy.H2PlsCool;
                
                fprintf( io, "\t%.4e", MAX2(0.,CoolHeavy.brems_cool_net) );
                debug_ctot += MAX2(0.,CoolHeavy.brems_cool_net);
                
                fprintf( io, "\t%.4e", CoolHeavy.heavfb );
                debug_ctot += CoolHeavy.heavfb;
                
                fprintf( io, "\t%.4e", CoolHeavy.eebrm );
                debug_ctot += CoolHeavy.eebrm;
                
                fprintf( io, "\t%.4e", CoolHeavy.tccool );
                debug_ctot +=  CoolHeavy.tccool;
                
                fprintf( io, "\t%.4e", CoolHeavy.cextxx );
                debug_ctot += CoolHeavy.cextxx;
                
                fprintf( io, "\t%.4e", CoolHeavy.expans );
                debug_ctot += CoolHeavy.expans;
                
                fprintf( io, "\t%.4e", CoolHeavy.cyntrn );
                debug_ctot += CoolHeavy.cyntrn;
                
                fprintf( io, "\t%.4e", CoolHeavy.colmet );
                debug_ctot += CoolHeavy.colmet;
                
                fprintf( io, "\t%.4e", thermal.dima );
                debug_ctot += thermal.dima;
                fprintf( io, " \n" );
                if( fabs( (debug_ctot - cool_total)/cool_total ) > 1e-10 )
                {
                        fprintf( ioQQQ , "PROBLEM with the SAVE EACH COOLING output\n" );
                        fprintf( ioQQQ , "PROBLEM One or more coolants have been lost, the sum of the reported cooling is %.4e\n", debug_ctot );
                        fprintf( ioQQQ , "PROBLEM The total cooling is %.4ee\n", cool_total );
                        fprintf( ioQQQ , "PROBLEM The difference is %.4e\n", cool_total - debug_ctot );
                        cdEXIT(EXIT_FAILURE);
                }
        }
        else if( strcmp(chJob,"COOL") == 0 )
                {
                /*>>chng 06 jun 06, change start of save to give same info as heating 
                 * as per comment by Yumihiko Tsuzuki */
                /* begin the print out with zone number, total heating and cooling */
                fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e", 
                                 radius.depth_mid_zone, 
                                 phycon.te, 
                                 heat_total, 
                                 cool_total );
 
                /* print only up to IPRINT, which is defined above */
                ip = MIN2( ip , IPRINT );
 
                /* now print the coolants 
                 * keep sign of coolant, for strong negative cooling 
                 * order is ion, wavelength, fraction of total */
                for( is=0; is < ip; is++ )
                {
                        if(is > 4 && (thermal.cooling[index[is]] < cset && thermal.heatnt[index[is]] < cset))
                                break;
                        fprintf( io, "\t%s %.1f\t%.7f", 
                                         thermal.chClntLab[index[is]], 
                                         thermal.collam[index[is]], 
                                         sign(csav[index[is]],sgnsav[index[is]]) );
                }
                fprintf( io, " \n" );
        }
        else
                TotalInsanity();
 
        /* finished, now free space */
        free(sgnsav);
        free(csav);
        free(index);
 
        return;
}