	program axion
c Agnes caught a missing T^4 on line with label 110
c 20 may 2004 - corrected T**4 on L line (label 110)
C Axion luminosity for C12 core accordind to Masayuki Nakagawa et al
C  1988, ApJ, 326, 241
	implicit double precision(a-h,o-z)
	real*8 pi,alphaa,logrho,logT,T,rho,Gamma,u,v,w,x,cosu,sinu,Fu1,Fu160,Fu5000,Fu180p,F,L,logFliquid,logFlattice,logFphonon
        pi=dacos(-1.0d0)
C ATTENTION:	alphaa=1e-26
 	alphaa=1.0d-26
C 1-9 are for 0.6 Msun C/O wd
	do 200 i=1,9
        if(i.eq.1)logT=8.1742
        if(i.eq.1)logrho=6.259
        if(i.eq.2)logT=7.8279
        if(i.eq.2)logrho=6.5050
        if(i.eq.3)logT=7.5792
        if(i.eq.3)logrho=6.5410
        if(i.eq.4)logT=7.0568
        if(i.eq.4)logrho=6.5569
        if(i.eq.5)logT=6.9198
        if(i.eq.5)logrho=6.5587
        if(i.eq.6)logT=6.6244
        if(i.eq.6)logrho=6.5611
        if(i.eq.7)logT=6.1452
        if(i.eq.7)logrho=6.5629
        if(i.eq.8)logT=5.9823
        if(i.eq.8)logrho=6.5631
        if(i.eq.9)logT=5.7627
        if(i.eq.9)logrho=6.5632
C 10-18 are for 1.0 Msun C/O wd
        if(i.eq.10)logT=8.0126
        if(i.eq.10)logrho=7.4730
        if(i.eq.11)logT=7.8034
        if(i.eq.11)logrho=7.5020
        if(i.eq.12)logT=7.2794
        if(i.eq.12)logrho=7.5192
        if(i.eq.13)logT=6.9477
        if(i.eq.13)logrho=7.5234
        if(i.eq.14)logT=6.8036
        if(i.eq.14)logrho=7.5246
        if(i.eq.15)logT=6.6585
        if(i.eq.15)logrho=7.5255
        if(i.eq.16)logT=6.3525
        if(i.eq.16)logrho=7.5263
        if(i.eq.17)logT=5.8726
        if(i.eq.17)logrho=7.5265
        if(i.eq.18)logT=5.4889
        if(i.eq.18)logrho=7.5266
	T=10.0d0**logT
	rho=10.0d0**logrho
        Gamma=3.577/(T/1.0d8)*(rho/1.0d6)**0.33333
	u=2.0d0*pi*(logrho+5.0d0)/18.0d0
	v=-0.3423d0+2.2053d0/Gamma**(0.3333)-2.1415d0/Gamma**(0.6666)+1.3474d0/Gamma
	w=0.1906d0-5.6458d0/Gamma**(0.3333)+32.461d0/Gamma**(0.6666)+142.6d0/Gamma
	x=12.118d0-197.1d0/Gamma**(0.333)+1253.6d0/Gamma**(0.6666)-2804.8d0/Gamma
	cosu=dcos(u)
	sinu=dsin(u)
C	write(*,*)u,v,w,x
	if (Gamma.gt.178d0) goto 50
	Fu1=2.7337d0/2.0d0-0.8648d0*cosu-0.2367d0*(cosu)**2-0.0715d0*(cosu)**3-0.0477d0*(cosu)**4
        Fu1=Fu1-0.345d0*sinu-0.0135d0*sinu**2+0.0132d0*sinu**3+0.0022d0*sinu**4-0.1395d0*u-1.3894d0
	Fu160=3.1029d0/2-1.0355*cosu-0.247*(cosu)**2-0.0551*(cosu)**3-0.0558*(cosu)**4
        Fu160=Fu160-0.3332*sinu+0.0271*sinu**2+0.005*sinu**3-0.0026*sinu**4-0.3146*u-1.3409d0
C	write(*,*)Fu1,Fu160
	logFliquid=v*Fu1+(1.0d0-v)*Fu160
	F=10.0d0**logFliquid
	goto 110
50	Fu180=3.565d0/2-1.1495*cosu-0.2968*(cosu)**2-0.0738*(cosu)**3-0.0694*(cosu)**4
        Fu180=Fu180-0.4141*sinu+0.0072*sinu**2+0.0007*sinu**3-0.003*sinu**4-0.3664*u-1.7456d0
	Fu5000=3.621d0/2-1.1746*cosu-0.297*(cosu)**2-0.0745*(cosu)**3-0.0696*(cosu)**4
	Fu5000=F15000-0.4223*sinu+0.012*sinu**2+0.0007*sinu**3-0.0033*sinu**4-0.4079*u-1.4663d0
	Fu180p=3.8289d0/2-1.1987*cosu-0.3269*(cosu)**2-0.0939*(cosu)**3-0.0787*(cosu)**4
	Fu180p=Fu180p-0.5103*sinu-0.0101*sinu**2-0.0034*sinu**3-0.0018*sinu**4-0.4501*u-1.9543d0
C	write(*,*)Fu180,Fu5000,Fu180p
	logFlattice=w*Fu180+(1.0d0-w)*Fu5000
	logFphonon=x*Fu180p
	F=10.0d0**logFlattice+10.0d0**logFphonon
c 20 may 2004 - corrected T**4 on next line
110	L=3.24d23*alphaa*(T/1.0d7)**4*rho*F
	write(*,100)Gamma,L
100	format(1x,f10.4,1x, 1pe12.6)
200	continue
	stop
	end
