/*CoolCalc compute calcium cooling */
#include "cddefines.h"
#include "coolheavy.h"
#include "ionfracs.h"
#include "ca2mly.h"
#include "taulines.h"
#include "tepowers.h"
#include "doppvel.h"
#include "phycon.h"
#include "dcala.h"
#include "caline.h"
#include "opac.h"
#include "edsqte.h"
#include "pmpcah.h"
#include "rfield.h"
#include "fivel.h"
#include "ligbar.h"
#include "makecs.h"
#include "coladd.h"
#include "level2.h"
#include "putcs.h"
#include "level3.h"
#include "pop3.h"
#include "coolmetals.h"
#include "pca2ex.h"

void CoolCalc()
{
	float p2 ;
	double  a21, 
	  a31, 
	  a41, 
	  a42, 
	  a51, 
	  a52, 
	  a53, 
	  c21, 
	  Ca2pop[5] ,
	  cs, 
	  cs2s2p, 
	  cs2s3p	, 
	  cs01, 
	  cs02, 
	  cs12, 
	  cs14, 
	  cs15, 
	  d41, 
	  d42, 
	  d51, 
	  d52, 
	  d53, 
	  hlgam, 
	  op41, 
	  op51, 
	  opckh, 
	  opcxyz, 
	  p21, 
	  p3, 
	  p31, 
	  p41, 
	  p51, 
	  r21, 
	  r31, 
	  r41, 
	  r42, 
	  r51, 
	  r52, 
	  r53;
	static double gCa2[5]={2.,4.,6.,2.,4.};
	static double exCa2[4]={13650.2,60.7,11480.6,222.9};
	static float opCax = 0.f;
	static float opCay = 0.f;
	static float opCaz = 0.f;
	
#	ifdef DEBUG_FUN
	fputs( "<+>CoolCalc()\n", debug_fp );
#	endif

	/* Ca I 4228 */
	MakeCS(&TauLines[ipCaI4228]);
	level2(&TauLines[ipCaI4228]);

	/* photoionization of evcited levels by Ly-alpha */
	hlgam = rfield.otslin[ HydroLines[0][IP2P][IP1S].ipCont -1];
	p51 = 1.7e-18*hlgam;
	p41 = 8.4e-19*hlgam;
	p31 = 7.0e-18*hlgam;
	p21 = 4.8e-18*hlgam;

	/* spontaneous decays
	 * frist two trans prob from 
	 * >>refer Zeippen, C.J. 1990, A&A, 229, 248 */
	a21 = 1.02*TauLines[ipT7324].Pesc;
	a31 = 1.05*TauLines[ipT7291].Pesc;
	a41 = 1.4e8*TauLines[ipT3969].Pesc;
	a51 = 1.4e8*TauLines[ipT3934].Pesc;
	a42 = 7.9e6*TauLines[ipT8662].Pesc;
	a52 = 8.2e5*TauLines[ipT8498].Pesc;
	a53 = 7.48e6*TauLines[ipT8542].Pesc;

	/* destruction of IR triplet by continuous opacities */
	opcxyz = opac.opac[ TauLines[ipT7324].ipCont -1];

	/* opcxyz = opac(icaxyz) */
	if( opcxyz > 0. )
	{
		d52 = 5.6*opcxyz/(opcxyz + opCax)*(1. - TauLines[ipT8498].Pesc);
		d53 = 5.6*opcxyz/(opcxyz + opCay)*(1. - TauLines[ipT8542].Pesc);
		d42 = 5.6*opcxyz/(opcxyz + opCaz)*(1. - TauLines[ipT8662].Pesc);
	}
	else
	{
		d52 = 0.;
		d53 = 0.;
		d42 = 0.;
	}

	/* near UV dest of KH by background continuum */
	opckh = opac.opac[ TauLines[ipT3969].ipCont -1];

	/* opckh = opac(icakh) */
	if( opckh > 0. )
	{
		op51 = xIonFracs[2][19]*3.89e-7/DoppVel.doppler[19];
		d51 = 5.6*opckh/(opckh + op51);
		op41 = xIonFracs[2][19]*1.96e-7/DoppVel.doppler[19];
		d41 = 5.6*opckh/(opckh + op41);
	}
	else
	{
		op51 = 0.;
		d51 = 0.;
		op41 = 0.;
		d41 = 0.;
	}
	/* net rates */
	r21 = p21 + a21;
	r31 = p31 + a31;
	r41 = a41 + p41 + d41;
	r51 = a51 + p51 + d51;
	r42 = a42 + d42;
	r52 = a52 + d52;
	r53 = a53 + d53;
	cs14 = 0.923*TePowers.te10*TePowers.te10;
	cs15 = cs14*2.;
	TauLines[ipT3969].cs = (float)cs14;
	TauLines[ipT3934].cs = (float)cs15;

	/* following used to correct rec contribution
	 * fcakh = a51 / ( a51 + eden*1.5e-5 / sqrte ) 
	 * cs 1-2 from 
	 * >>refer Saraph, H.E. 1970, J.Phys. B, 3, 952
	 * other 
	 * >>refer Chidichimo, M.C. 1981, J.Phys. B, 14, 4149 */
	fivel(gCa2,exCa2,5.8,8.6,cs14,cs15,20.6,22.9,9.8,3.4,44.4,1.0,
	  r21,r31,r41,r51,0.,r42,r52,0.,r53,0.,Ca2pop,xIonFracs[2][19]);

	/* CDSQTE = 8.629E-6*EDEN/SQRTE */
	c21 = 5.8/4.*edsqte.cdsqte;

	/* remember largest ratio of Ly-al removal to total */
	if( xIonFracs[2][19] > 0. )
	{
		ca2mly.Ca2RmLya = (float)MAX2(ca2mly.Ca2RmLya,p21/(p21+a21+c21));
	}
	caline.Cak = (float)(Ca2pop[4]*a51*5.06e-12);
	caline.Cah = (float)(Ca2pop[3]*a41*5.01e-12);
	caline.Cax = (float)(Ca2pop[4]*a52*2.34e-12);
	caline.Cay = (float)(Ca2pop[4]*a53*2.33e-12);
	caline.Caz = (float)(Ca2pop[3]*a42*2.30e-12);
	caline.Caf1 = (float)(Ca2pop[2]*a31*2.73e-12);
	caline.Caf2 = (float)(Ca2pop[1]*a21*2.72e-12);
	pca2ex.popca2ex = (float)(Ca2pop[1] + Ca2pop[2] + Ca2pop[3] + Ca2pop[4]);

	/* this is the total cooling due to the model atom */
	CoolHeavy.Cair = caline.Cax + caline.Cay + caline.Caz;
	CoolHeavy.c7306 = caline.Caf1 + caline.Caf2;
	CoolHeavy.Cakh = caline.Cak + caline.Cah;

	/* add the CaII lines to the cooling stack */
	coladd("Ca 2",7306,CoolHeavy.c7306);
	coladd("Ca 2",8400,CoolHeavy.Cair);
	coladd("Ca 2",3800,CoolHeavy.Cakh);

	/* level populations that will be used for excited state photoionization */
	dcala.dstCala = (float)(Ca2pop[4]*p51 + Ca2pop[3]*p41);
	dcala.dCakh = (float)(dcala.dstCala*5.03e-12);
	dcala.dCaf12 = (float)((Ca2pop[2]*p31 + Ca2pop[1]*p21)*2.31e-12);
	opCax = (float)(Ca2pop[1]*1.13e-8/DoppVel.doppler[19]);
	opCay = (float)(Ca2pop[2]*6.87e-8/DoppVel.doppler[19]);
	opCaz = (float)(Ca2pop[1]*5.74e-8/DoppVel.doppler[19]);

	/* total rate Lalpha destroys CaII,
	 * this is only used in ioncali to increase ionization rate by
	 * adding it to the ct vector */
	if( xIonFracs[2][19] > 0. )
	{
		dcala.dstCala = (float)(
			(dcala.dstCala + dcala.dCaf12/2.31e-12)/xIonFracs[2][19]);
		/* take average of old and new */
		dcala.dstCala = (float)((dcala.dstCala + dcala.oldcla)/2.);
		dcala.oldcla = dcala.dstCala;
		{

			enum {DEBUG=FALSE};
			if( DEBUG )
			{
				fprintf(ioQQQ," hlgam is %e\n", hlgam);
			}
		}
	}
	else
	{
		dcala.dstCala = 0.;
	}
	dcala.Ca3d = (float)(Ca2pop[1] + Ca2pop[2]);
	dcala.Ca4p = (float)(Ca2pop[3] + Ca2pop[4]);

	/* incl stimulated emission for Calcium II 5-level atom */
	TauLines[ipT3934].PopOpc = (float)(Ca2pop[0] - Ca2pop[4]/2.);
	TauLines[ipT3934].PopHi = (float)Ca2pop[4];
	TauLines[ipT3934].PopLo = (float)Ca2pop[0];
	TauLines[ipT3969].PopOpc = (float)(Ca2pop[0] - Ca2pop[3]);
	TauLines[ipT3969].PopHi = (float)Ca2pop[3];
	TauLines[ipT3969].PopLo = (float)Ca2pop[0];

	pmpcah.tpcah[0] = TauLines[ipT3969].TauIn;
	TauLines[ipT8498].PopOpc = (float)(Ca2pop[1] - Ca2pop[4]);
	TauLines[ipT8498].PopHi = (float)Ca2pop[4];
	TauLines[ipT8498].PopLo = (float)Ca2pop[1];
	TauLines[ipT8542].PopOpc = (float)(Ca2pop[2] - Ca2pop[4]*1.5);
	TauLines[ipT8542].PopHi = (float)Ca2pop[4];
	TauLines[ipT8542].PopLo = (float)Ca2pop[2];
	TauLines[ipT8662].PopOpc = (float)(Ca2pop[1] - Ca2pop[3]*2.);
	TauLines[ipT8662].PopHi = (float)Ca2pop[3];
	TauLines[ipT8662].PopLo = (float)Ca2pop[1];
	TauLines[ipT7291].PopOpc = xIonFracs[2][19];
	TauLines[ipT7291].PopHi = 0.;
	TauLines[ipT7291].PopLo = xIonFracs[2][19];
	TauLines[ipT7324].PopOpc = xIonFracs[2][19];
	TauLines[ipT7324].PopHi = 0.;
	TauLines[ipT7324].PopLo = xIonFracs[2][19];

	/* Ca IV 3.2 micron; data from
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * Y(ik) from 
	 * >>refer Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209 */
	if( phycon.te <= 1e5 )
	{
		cs = MAX2(1.0,8.854e-3*TePowers.sqrte);
	}
	else if( phycon.te < 2.512e5 )
	{
		cs = 2.8;
	}
	else
	{
		cs = 641.1/(TePowers.te30*TePowers.te10*TePowers.te02*TePowers.te02/
		  TePowers.te003);
	}
	PutCS(cs,&TauLines[ipTCa3]);
	level2(&TauLines[ipTCa3]);

	/* [Ca V] IR 4.16, 11.47 micron; A from
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * cs from 
	 * >>refer Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
	 * >>chng 96 july 16, big changes in cs */
	cs = MIN2(3.3,0.392*TePowers.te20/TePowers.te005/TePowers.te003);
	cs = MAX2(2.2,cs);
	PutCS(cs,&TauLines[ipTCa4]);

	/* >>chng 96 aug 02, following had error in te dep in ver 90.01 */
	cs = MIN2(0.93,0.162*TePowers.te10*TePowers.te05*TePowers.te003*
	  TePowers.te001);
	cs = MAX2(0.67,cs);
	PutCS(cs,&TauLines[ipTCa12]);

	cs = MIN2(0.97,0.0894*TePowers.te20*TePowers.te01*TePowers.te005);
	cs = MAX2(0.60,cs);
	PutCS(cs,&TauDummy);

	level3(&TauLines[ipTCa4],&TauLines[ipTCa12],&TauDummy);

	/* Ca V lines from 1d, 1s; A from 
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * cs from 
	 * >>refer Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
	 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
	cs01 = MIN2(4.1,0.533*TePowers.te20/TePowers.te01);
	cs01 = MAX2(2.8,cs01);
	cs02 = MIN2(0.87,5.22e-03*TePowers.sqrte);
	p3 = pop3(9.,5.,1.,cs01,cs02,1.35,2.326,23.2,3.73,2.57e4,3.60e4,
	  &p2,xIonFracs[5][19],0.);

	CoolHeavy.c3997 = p3*3.73*4.98e-12;
	CoolHeavy.c2414 = p3*23.1*8.245e-12;
	CoolHeavy.Ca6087 = p2*0.426*3.268e-12;
	CoolHeavy.c5311 = p2*1.90*3.747e-12;

	coladd("Ca 5",3997,CoolHeavy.c3997);
	coladd("Ca 5",2414,CoolHeavy.c2414);
	coladd("Ca 5",6087,CoolHeavy.Ca6087);
	coladd("Ca 5",5311,CoolHeavy.c5311);

	/* Ca VII lines from 1d, 1s
	 * all cs from 
	 * >>refer Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
	 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
	cs01 = MIN2(4.4,22.25/(TePowers.te20/TePowers.te02/TePowers.te02));
	cs01 = MAX2(3.5,cs01);

	cs12 = MIN2(1.20,0.303*TePowers.te30*TePowers.te03);
	cs12 = MAX2(0.62,cs12);

	cs02 = MIN2(0.959,7.889/(TePowers.te20*TePowers.te05/TePowers.te01));
	cs02 = MAX2(0.50,cs02);

	p3 = pop3(9.,5.,1.,cs01,cs02,cs12,3.124,30.4,6.81,2.91e4,3.90e4,
	  &p2,xIonFracs[7][19],0.);

	CoolHeavy.Ca3688 = p3*6.81*5.40e-12;
	CoolHeavy.Ca2112 = p3*30.4*9.42e-12;
	CoolHeavy.Ca5620 = p2*2.15*3.548e-12;
	CoolHeavy.Ca4941 = p2*0.974*4.037e-12;
	coladd("Ca 7",3688,CoolHeavy.Ca3688);
	coladd("Ca 7",2112,CoolHeavy.Ca2112);
	coladd("Ca 7",5620,CoolHeavy.Ca5620);
	coladd("Ca 7",4941,CoolHeavy.Ca4941);

	/* all cs from 
	 * >>refer Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
	 * [Ca VII] 4.09, 6.15 mic 3P lines */
	cs = MIN2(5.354,0.406*TePowers.te20*TePowers.te03*TePowers.te01);
	cs = MAX2(3.702,cs);
	PutCS(cs,&TauLines[ipCa0741]);

	cs = MIN2(1.59,0.183*TePowers.te20);
	cs = MAX2(1.153,cs);
	PutCS(cs,&TauLines[ipCa0761]);

	cs = MIN2(1.497,0.0917*TePowers.te20*TePowers.te05* TePowers.te01);
	cs = MAX2(1.005,cs);
	PutCS(cs,&TauDummy);

	/* level3(  t10,t21,t20) */
	level3(&TauLines[ipCa0761],&TauLines[ipCa0741],&TauDummy);

	/* [Ca VIII]  2.32 microns, cs 
	 * >>refer Saraph, H.E., & Storey, P.J. A&AS, 115, 151 */
	cs = MIN2(6.75,22.04/(TePowers.te10*TePowers.te02*TePowers.te005));
	PutCS(cs,&TauLines[ipCa08232]);
	level2(&TauLines[ipCa08232]);

	/* [Ca 12] 3328.78A, cs from 
	 * >>refer Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
	cs = MIN2(0.172,0.0118*TePowers.te20*TePowers.te01);
	cs = MAX2(0.10,cs);
	PutCS(cs,&TauLines[ipCa12333]);
	level2(&TauLines[ipCa12333]);

	/* Li seq 2s2p 2s3p, 2s2p as two separate lines
	 * >>refer Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(20,&TauLines[ipTCa302],&TauLines[ipTCa19],&cs2s2p,&cs2s3p);

	PutCS(cs2s2p,&TauLines[ipTCa302]);
	level2(&TauLines[ipTCa302]);

	/* funny factor (should have been 0.5) due to energy change */
	PutCS(cs2s2p*0.439,&TauLines[ipTCa345]);
	level2(&TauLines[ipTCa345]);

	PutCS(cs2s3p,&TauLines[ipTCa19]);
	level2(&TauLines[ipTCa19]);


#	ifdef DEBUG_FUN
	fputs( " <->CoolCalc()\n", debug_fp );
#	endif
	return;
}