/*CoolSulf compute sulphur cooling */
#include "cddefines.h"
#include "coolheavy.h"
#include "sii.h"
#include "taulines.h"
#include "tepowers.h"
#include "phycon.h"
#include "embesq.h"
#include "poplevls.h"
#include "tcool.h"
#include "logte.h"
#include "ionfracs.h"
#include "s2cs.h"
#include "ligbar.h"
#include "coladd.h"
#include "level2.h"
#include "putcs.h"
#include "beseq.h"
#include "level3.h"
#include "pop3.h"
#include "popexc.h"
#include "coolmetals.h"
#include "expion.h"
#include "fivel.h"

void CoolSulf()
{
	double cs, 
	  cs2s2p, 
	  cs2s3p, 
	  cs01, 
	  cs02, 
	  cs12, 
	  cs13, 
	  cs14, 
	  cs15, 
	  cs23, 
	  cs24, 
	  cs25, 
	  cs34, 
	  cs35, 
	  cs45, 
	  p[5];
	float 
	  pop2;
	static double gS2[5]={4.,4.,6.,2.,4.};
	static double exS2[4]={14851.9,31.5,9640.8,48.6};
	

#	ifdef DEBUG_FUN
	fputs( "<+>CoolSulf()\n", debug_fp );
#	endif

	/* Sulphur II [S II] 6731 + 6716, A 
	 * >>refer Mendoza, C., & Zeippen, C.J., 1982, MNRAS, 198, 127
	 * collision strength from 
	 *>>refer Ramsbottom, C.A., Bell, K.L., Stafford, R.P. 1996,
	 *>>refer At. Data Nucl. Data Tables, 63, 57 */
	S2cs(logte.alogte,&cs12,&cs13,&cs14,&cs15,&cs23,&cs24,&cs25,&cs34,
	  &cs35,&cs45);
	/* cs13 = MIN(4.95 , 3.2640828 + 1.9369099*sexp(te/30337.498) )
	 * cs12 = cs13 / 1.5
	 * cs15 = MIN( 1.59 , 0.10*te10*te10*te10/te01/te003/te001/te001)
	 * cs14 = cs15 * 0.5
	 * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45,
	 *  A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */
	fivel(gS2,exS2,cs12,cs13,cs14,cs15,cs23,cs24,cs25,cs34,cs35,cs45,
	  8.82e-4,2.60e-4,9.06e-2,0.225,3.35e-7,0.163,0.133,0.078,0.179,
	  1.03e-6,p,xIonFracs[2][15]);
	sii.S6733 = (float)(p[1]*8.82e-4*2.96e-12);
	sii.S6718 = (float)(p[2]*2.60e-4*2.962e-12);
	sii.S4070 = (float)(p[4]*0.225*4.89e-12);
	sii.S4078 = (float)(p[3]*0.0906*4.88e-12);
	sii.S10323 = (float)(p[4]*0.179*1.93e-12);
	sii.S10289 = (float)(p[4]*0.133*1.93e-12);
	sii.S10373 = (float)(p[3]*0.078*1.92e-12);
	sii.S10339 = (float)(p[3]*0.163*1.92e-12);
	CoolHeavy.c6731 = sii.S6733 + sii.S6718;
	CoolHeavy.c10330 = sii.S10323 + sii.S10289 + sii.S10373 + sii.S10339;
	tcool.dCooldT += CoolHeavy.c6731*(2.14e4*tcool.tsq1 - tcool.halfte);
	tcool.dCooldT += CoolHeavy.c10330*2.58*(3.53e4*tcool.tsq1 - tcool.halfte);
	coladd("S  2",6731,CoolHeavy.c6731);
	coladd("S  2",10330,CoolHeavy.c10330);
	coladd("S  2",4070,CoolHeavy.c10330*1.58);

	/* SII 1256 */
	cs = MIN2(8.46,-4.9416304+47.01064/logte.alogte);
	cs = MAX2(4.466,cs);
	PutCS(cs,&TauLines[ipT1256]);
	level2(&TauLines[ipT1256]);

	/* S III 18.7M, 33.6M, A 
	 * >>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 99 dec 22, to Tayal, S.S., and Gupta, G.P. 1999 ApJ 526, 544 */
	cs = MIN2(2.331,7.935*phycon.te/(TePowers.te10*TePowers.te03*TePowers.te003));
	PutCS(cs,&TauLines[ipTS34]);
	if( phycon.te <= 39811. )
	{
		cs = MIN2(5.78,3.114*TePowers.te03*TePowers.te03);
	}
	else
	{
		cs = 24.93/(TePowers.te10*TePowers.te03*TePowers.te01/TePowers.te001/
		  TePowers.te001);
	}
	PutCS(cs,&TauLines[ipTS19]);
	cs = MIN2(1.413,0.221*phycon.te*TePowers.te20/TePowers.te03*TePowers.te005);
	PutCS(cs,&TauDummy);
	level3(&TauLines[ipTS34],&TauLines[ipTS19],&TauDummy);

	/* S III O III-like lines, A from 
	 * >>refer Mendoza, C., & Zeippen, C.J. 1982, MNRAS, 199, 1025
	 * 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) */
	cs = MIN2(2.05,0.0821*TePowers.te30);
	CoolHeavy.c6312 = pop3(9.,5.,1.,7.98,1.14,cs,7.97e-2,0.807,2.22,
	  1.55e4,2.28e4,&pop2,xIonFracs[3][15],0.)*2.22*3.15e-12;
	/* folowing is 9532 + 9069 together (OIII-like) */
	CoolHeavy.c9532 = pop2*7.93e-2*2.11e-12;
	coladd("S  3",6312,CoolHeavy.c6312);
	coladd("S  3",9532,CoolHeavy.c9532);
	coladd("S  3",3722,CoolHeavy.c6312*0.59);

	/* S III 1194, data from 
	 * >>refer Ho, Y.K., & Henry,R.J.W. 1984, ApJ, 282, 816
	 * >>chng 97 may 17, to, about 2x larger than above
	 * >>refer Tayal, S.S. 1997, ApJ 481, 550 */
	if( phycon.te <= 3e4 )
	{
		cs = 12.04/(TePowers.te02*TePowers.te02);
	}
	else if( phycon.te > 3e4 && phycon.te <= 4e4 )
	{
		cs = 7.97;
	}
	else
	{
		cs = 55.42/(TePowers.te20/TePowers.te02*TePowers.te003);
	}
	PutCS(cs,&TauLines[ipT1194]);
	level2(&TauLines[ipT1194]);

	/* S III] 1713.12, 1728.94, cs from 
	 * >>refer Hayes, M.A., 1986, J Phys B 19, 1853.
	 * cs = MIN( 4.0 , 7.794 / (te10/te02/te001/te001)  )
	 * >>chng 97 may 17, about 20% smaller than before
	 * >>refer Tayal, S.S. 1997, ApJ 481, 550 */
	if( phycon.te <= 3e4 )
	{
		cs = 1.786*TePowers.te05*TePowers.te01*TePowers.te001;
	}
	else
	{
		cs = 9.392/TePowers.te10;
	}
	PutCS(cs,&TauLines[ipTS1720]);
	level2(&TauLines[ipTS1720]);

	/* S IV 10.5MI, CS and A FROM 
	 * >>refer Johnson, C.T., Kingston, A.E., Dufton, P.L. 1986, MNRAS, 220, 155
	 * >>chng 97 feb 14, error in cs below t = 10,000K
	 * >>chng 96 dec 19, to CS from 
	 * >>refer Saraph, H.E., Storey, P.J., & Tully, J.A. 1995, 5th International
	 * >>refer Colloquium on Atomic Spectra and Oscillator Strengths, ed. by
	 * >>refer W.-U L.  Tchang-Brillet, J.-F. Wyart, C.J. Zeippen,
	 * >>refer (Meudon: Publications de l'Observaroire de Paris), p.110
	 * above said to be A&A in press */
	if( phycon.te < 1e4 )
	{
		cs = 3.71*TePowers.te10/TePowers.te01;
	}
	else
	{
		cs = MIN2(8.5,19.472/(TePowers.te10/TePowers.te01));
	}
	PutCS(cs,&TauLines[ipTS11]);
	level2(&TauLines[ipTS11]);

	/* S IV 1406, cs from 
	 * >>refer Dufton, P.L., Hibbert, A., Kingston, A.E., & Doschek, G.A. 1982,
	 * >>refer ApJ, 257, 338 */
	cs = 16.67/(TePowers.te10*TePowers.te02);
	PutCS(cs,&TauLines[ipT1406]);
	level2(&TauLines[ipT1406]);

	/* S V Be-seq line, A=
	 * >>refer Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
	 * >>refer ed by D.R. Flower, (D. Reidel: Holland), 143
	 * all cs 
	 * >>refer Dufton, P.L., Hibbert, A., Keenan, F.P, Kingston, A.E., &
	 * >>refer Doschek, G.A. 1986, ApJ, 300, 448
	 * and 
	 * >>refer Dufton, P.L., & Kingston, A.E. 1984, J.Phys. B, 17, 3321 */
	cs = MIN2(1.58,35.372/(TePowers.te10*TePowers.te10*TePowers.te10));
	PutCS(cs,&TauLines[ipT1198]);
	/* cs01 = MIN(1.98, 29.625/(T**0.25)) */
	cs01 = MIN2(1.98,29.625/(TePowers.te20*TePowers.te05));
	/* cs02 = MIN(2.26, 13.477/(T**0.165)) */
	cs02 = MIN2(2.26,13.477/(TePowers.te10*TePowers.te03*TePowers.te03*
	  TePowers.te005));
	/* cs12 = MIN(7.59, 63.994/(T**0.197)) */
	cs12 = MIN2(7.59,63.994/(TePowers.te20/TePowers.te003));
	/* beseq(cs23,cs24,cs34,tarray,a41)
	 * call beseq( 1.98 , 2.26 , 7.59, t1198 ,.066) */
	beseq(cs01,cs02,cs12,&TauLines[ipT1198],.066);
	embesq.em1198 = (float)(PopLevls.PopLevels[3]*0.066*1.66e-11);

	/* S V 786
	 * cs from 
	 * >>refer Dufton, P.L., & Kingston, A.E. 1984, J.Phys. B, 17, 3321 */
	PutCS(8.3,&TauLines[ipT786]);
	level2(&TauLines[ipT786]);

	/* S VI, Li sequence, 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 Dufton, P.L., & Kingston, A.E. 1987, J.Phys. B, 20, 3899 */
	cs = 11.9;
	PutCS(cs*0.667,&TauLines[ipT933]);
	PutCS(cs*0.333,&TauLines[ipT944]);
	PutCS(1.0,&TauDummy);

	level3(&TauLines[ipT944],&TauDummy,&TauLines[ipT933]);

	/* S VIII 9913
	 * >>refer Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
	cs = MIN2(0.291,0.0289*TePowers.te20/TePowers.te01* TePowers.te001);
	cs = MAX2(0.192,cs);
	PutCS(cs,&TauLines[ipfsS810]);

	level2(&TauLines[ipfsS810]);

	/* S IX 1715, 1987, 
	 * >>refer Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */
	CoolHeavy.c1715 = popexc(0.4,9.,5.,80.,8.39e4,xIonFracs[9][15])*
	  1.16e-11;
	coladd("S  9",1715,CoolHeavy.c1715);

	/* S IX 1.25, 3.78 microns
	 * cs 
	 * >>refer Butler, K., & Zeippen, C.J. 1994, A&AS, 108, 1 */
	cs = MIN2(3.10,277.2/(TePowers.sqrte/TePowers.te003));
	PutCS(cs,&TauLines[ipfsS912]);

	cs = MIN2(.50,7.940/(TePowers.te10*TePowers.te10*TePowers.te10));
	PutCS(cs,&TauLines[ipfsS938]);

	cs = MIN2(1.1,195.72/(TePowers.sqrte*TePowers.te03*TePowers.te03*
	  TePowers.te01*TePowers.te003*TePowers.te001));
	PutCS(cs,&TauDummy);

	level3(&TauLines[ipfsS912],&TauLines[ipfsS938],&TauDummy);

	/* S X 1213, 1197, 
	 * >>refer Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */
	CoolHeavy.c1213 = popexc(0.28,4.,6.,6.,1.19e5,xIonFracs[10][15])*
	  1.64e-11;
	coladd("S 10",1213,CoolHeavy.c1213);

	/* S XI 1.91, 1.40 microns, cs 
	 * >>refer Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
	cs = MIN2(0.22,0.0151*TePowers.te20*TePowers.te03*TePowers.te001*
	  TePowers.te001);
	PutCS(cs,&TauLines[ipfsS1119]);
	cs = MIN2(0.6,0.0168*TePowers.te10*TePowers.te10*TePowers.te10*
	  TePowers.te01);
	PutCS(cs,&TauLines[ipfsS1114]);
	cs = MIN2(0.15,1.5e-3*TePowers.te30*TePowers.te10*TePowers.te005);
	PutCS(cs,&TauDummy);
	level3(&TauLines[ipfsS1119],&TauLines[ipfsS1114],&TauDummy);

	/* S XI 1615, 1826, 
	 * >>refer Kafatos, M., & Lynch, J.P. 1980, ApJS, 42, 611 */
	CoolHeavy.c1826 = popexc(0.3,9.,5.,100.,8.92e4,xIonFracs[11][15])*
	  1.23e-11;
	coladd("S 11",1615,CoolHeavy.c1826);

	/* S XII 7611.2, cs from 
	 * >>refer Oliva, E., Pasquali, A., & Reconditi, M. 1996, A&A, 305, 210
	 * >>refer Froese Fischer, C. 1983, J.Phys. B, 16, 157
	 * >>refer Saha, H.P., & Trefftz, E. 1983, SoPh, 87, 233 */
	PutCS(0.135,&TauLines[ipfsS1207]);
	level2(&TauLines[ipfsS1207]);

	/* S XII 520, actually group of four intercombination lines
	 * all atomic data from 
	 * >>refer Saha, H.P., & Trefftz, E. 1983, SoPh, 87, 233 */
	CoolHeavy.c520 = 0.03*1.65e-16*expion(2.8e5,xIonFracs[12][15]);
	coladd("S 12",520,CoolHeavy.c520);

	/* S XIII 488.4, 1909-like */
	CoolHeavy.c488 = popexc(0.11,1.,9.,3.3e5,2.948e5,xIonFracs[13][15])*
	  4.07e-11;
	coladd("S 13",488,CoolHeavy.c488);

	/* li seq 2s2p and 2s3p, S XIV 418, 446
	 * >>refer Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
	ligbar(16,&TauLines[ipTSu418],&TauLines[ipTSu30],&cs2s2p,&cs2s3p);
	PutCS(cs2s2p,&TauLines[ipTSu418]);
	PutCS(cs2s2p*0.5,&TauLines[ipTSu446]);
	PutCS(1.0,&TauDummy);
	level3(&TauLines[ipTSu446],&TauDummy,&TauLines[ipTSu418]);

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


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