cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*CoolPunch punch coolants */ 00004 #include "cddefines.h" 00005 #include "thermal.h" 00006 #include "dynamics.h" 00007 #include "radius.h" 00008 #include "conv.h" 00009 #include "phycon.h" 00010 #include "punch.h" 00011 00012 /* this is limit to number of coolants to print out */ 00013 static const int IPRINT = 100; 00014 00015 /*CoolPunch punch coolants */ 00016 void CoolPunch(FILE * io) 00017 { 00018 long int i, 00019 ip, 00020 is; 00021 00022 int nFail; 00023 00024 double cset, 00025 cool_total, 00026 heat_total; 00027 00028 realnum 00029 *csav, 00030 *sgnsav; 00031 long int *index; 00032 00033 DEBUG_ENTRY( "CoolPunch()" ); 00034 00035 /* cannot do one-time init since thermal.ncltot can change */ 00036 index = (long int *)CALLOC((size_t)thermal.ncltot,sizeof(long int)); 00037 csav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum)); 00038 sgnsav = (realnum *)CALLOC((size_t)thermal.ncltot,sizeof(realnum)); 00039 00040 cool_total = thermal.ctot; 00041 heat_total = thermal.htot; 00042 00043 /* >>chng 06 mar 17, comment out following block and replace with this 00044 * removing dynamics heating & cooling and report only physical 00045 * heating and cooling 00046 * NB the heating and cooling as punched no longer need be 00047 * equal for a converged model */ 00048 cool_total -= dynamics.Cool; 00049 heat_total -= dynamics.Heat; 00050 # if 0 00051 if(dynamics.Cool > dynamics.Heat) 00052 { 00053 cool_total -= dynamics.Heat; 00054 heat_total -= dynamics.Heat; 00055 } 00056 else 00057 { 00058 cool_total -= dynamics.Cool; 00059 heat_total -= dynamics.Cool; 00060 } 00061 # endif 00062 00063 /* cset will be weakest cooling to consider 00064 * WeakHeatCool set with 'set weakheatcool' command 00065 * default is 0.05 */ 00066 cset = cool_total*punch.WeakHeatCool; 00067 00068 /* first find all strong lines, both + and - sign */ 00069 ip = thermal.ncltot; 00070 00071 for( i=0; i < ip; i++ ) 00072 { 00073 csav[i] = (realnum)(MAX2(thermal.cooling[i],thermal.heatnt[i])/ 00074 cool_total); 00075 00076 /* save sign to remember if heating or cooling line */ 00077 if( thermal.heatnt[i] == 0. ) 00078 { 00079 sgnsav[i] = 1.; 00080 } 00081 else 00082 { 00083 sgnsav[i] = -1.; 00084 } 00085 } 00086 00087 /* order strongest to weakest */ 00088 /* now sort by decreasing importance */ 00089 /*spsort netlib routine to sort array returning sorted indices */ 00090 spsort( 00091 /* input array to be sorted */ 00092 csav, 00093 /* number of values in x */ 00094 ip, 00095 /* permutation output array */ 00096 index, 00097 /* flag saying what to do - 1 sorts into increasing order, not changing 00098 * the original routine */ 00099 -1, 00100 /* error condition, should be 0 */ 00101 &nFail); 00102 00103 /* warn if tcovergence failure occurred */ 00104 if( !conv.lgConvTemp ) 00105 { 00106 fprintf( io, "#>>>> Temperature not converged.\n" ); 00107 } 00108 else if( !conv.lgConvEden ) 00109 { 00110 fprintf( io, "#>>>> Electron density not converged.\n" ); 00111 } 00112 else if( !conv.lgConvIoniz ) 00113 { 00114 fprintf( io, "#>>>> Ionization not converged.\n" ); 00115 } 00116 else if( !conv.lgConvPres ) 00117 { 00118 fprintf( io, "#>>>> Pressure not converged.\n" ); 00119 } 00120 00121 /*>>chng 06 jun 06, change start of punch to give same info as heating 00122 * as per comment by Yumihiko Tsuzuki */ 00123 /* begin the print out with zone number, total heating and cooling */ 00124 fprintf( io, "%.5e\t%.4e\t%.4e\t%.4e", 00125 radius.depth_mid_zone, 00126 phycon.te, 00127 heat_total, 00128 cool_total ); 00129 00130 /* print only up to IPRINT, which is defined above */ 00131 ip = MIN2( ip , IPRINT ); 00132 00133 /* now print the coolants 00134 * keep sign of coolant, for strong negative cooling 00135 * order is ion, wavelength, fraction of total */ 00136 for( is=0; is < ip; is++ ) 00137 { 00138 if(is > 4 && (thermal.cooling[index[is]] < cset && thermal.heatnt[index[is]] < cset)) 00139 break; 00140 fprintf( io, "\t%s %.1f\t%.7f", 00141 thermal.chClntLab[index[is]], 00142 thermal.collam[index[is]], 00143 sign(csav[index[is]],sgnsav[index[is]]) ); 00144 } 00145 fprintf( io, " \n" ); 00146 00147 /* finished, now free space */ 00148 free(sgnsav); 00149 free(csav); 00150 free(index); 00151 return; 00152 }