c  Program to calculate  spectrum of dust emission
c  from a circumstellar envelope and to convolve with
c  IRAS filters
c
c  Simple model - use Draine-Lee model for grains
c
c  The input file, called 'input', contains the following on one
c  line:
c  effective temperature of star, K
c  stellar luminosity, solar luminosities
c  outer radius of dust shell, cm
c  inner radius of dust shell, cm
c  dust loss loss rate in solar masses per year
c  outflow speed in km/sec
c  distance in parsecs
c  a flag for silicate or graphite grains
c  the grain radius in cm
c  the grain material density
c  the grain emissivity index at long wavelengths
c  the fraction of the mass in large grains
c  the radius of the large grains
c
	implicit double precision (a-h,o-z)
c r are the radial zones
	dimension flux(136),freq(136),fh(136),fjy(136),xint(136)
	dimension xlam(136),q(136),qsi(136),qc(136),filter(4,136)
	dimension f0(4),fluxtrue(4),fluxobs(4),cor(4)
	open(1,file='output', status='old')
	rewind 1
	open(2,file='input',status='old')
	rewind 2
	open(3,file='spect.dat',status='old')
	rewind 3
	open(4,file='star.dat',status='old')
	rewind 4
	open(8,file='irasmod.dat',status='old')
	rewind 8
	open(10,file='qfin.dat',status='old')
	open(11,file='filter.dat',status='old')
	open(12,file='flux.dat',status='old')
	open(13,file='color.dat',status='old')
	rewind 10
	rewind 11
	rewind 12
	rewind 13
	pi=3.1415926535
	sigma=5.66956d-5
	solarl=3.826d33
	solarm=1.989d33
	secyr=3.1558d7
	h=6.6262d-27
	c=2.9979248d10
	xk=1.38062d-16
	ag=2.d-5
	parsec=3.086d18
	xjansky=1.d23
	flxcon=xjansky/4./pi/parsec**2
	angcon=180.*3600./pi/parsec
	f0(1)=c/12.d-4
	f0(2)=c/25.d-4
	f0(3)=c/60.d-4
	f0(4)=c/1.d-2
	np=136
c  
	do 55 i=1,np
	read(11,*) dum,filter(1,i),filter(2,i),filter(3,i),filter(4,i)
55	continue
	read(2,*) tstar,starl,rmax,rmin,xml,v0,dist,igrtype,
     1  ag,rhog,beta,f,agl
	chi=(xml/1.d-7)*(15./v0)*(1.d4/starl)
1002	continue
c
c  Calculate radiation pressure optical depth
c	taum=1.d7*solarm*xml*v0*c/secyr/starl/solarl
	taum=1.d7*solarm*xml/secyr/starl/solarl
	taum=taum*v0*c
	write(1,217) taum
217	format('ratio of wind and photon momenta = ',1x,1pe10.3)
	write(1,5)
c
	xmll=f*xml
	xml=(1.-f)*xml
	i=0
41	i=i+1
  	read(10,*,end=40) ndum,xlam(i),qsi(i),qc(i)
	go to 41
40	continue	
	if(igrtype.eq.0) then
	tcond=1300.
	else
	tcond=1500.
	end if
	do 44 i=1,np
	if(igrtype.eq.0) then
	q(i)=qsi(i)
	else
	q(i)=qc(i)
	end if
44	continue
	xmg=4.*pi*ag**3*rhog/3.
	xmgl=4.*pi*agl**3*rhog/3.
	rstar=starl*solarl/4./pi/tstar**4/sigma
	rstar=dsqrt(rstar)
	flxcon=flxcon/dist**2
	angcon=angcon/dist
c  Calculate inner radius; grain temperature = 1300 K (silicates)
c  1500 K (graphite)
c
	if(rmin.gt.0.) go to 321
 	rmin=0.5*rstar*(tstar/tcond)**2.5
	if(rstar.ge.rmin) then
	print *,'stellar radius larger than inner radius of dust shell'
	print *, 'inner radius = ',rmin
	print *, 'stellar radius = ',rstar
	go to 1000
	end if
321	continue
	write(1,1) rstar,tstar,starl
1	format('stellar radius = ',1pd10.3,/,2x, 'stellar
     1  temperature = ',0pf7.1,/,2x,'stellar luminosity = ',
     2  1pe10.3, 'solar luminosities')
	write(1,11) rmax,rmin
11	format('maximum radius = ',1pe10.3,' minimum = ',1pe10.3)
	xnd0=4.993d19*xml/xmg/v0
	xnd0l=4.993d19*xmll/xmgl/v0
	write(1,21) xnd0,xnd0l
21	format('central dust densities',2(1x,1pe10.3))
	write(1,5)
	if(igrtype.eq.0) then
	write(1,22) ag,rhog
22	format('silicate grains: radius = ',1pe10.3,' cm',
     1  'density = ',0pf4.2,' gm/cm**3')
	write(1,5)
	end if
	if(igrtype.eq.1) then
	write(1,23) ag,rhog
23	format('graphite grains: radius = ',1pe10.3,' cm',
     1  'density = ',0pf4.2,' gm/cm**3')
	write(1,5)
	end if
	arstar=rstar*angcon
	armin=rmin*angcon
	armax=rmax*angcon
	write(1,18) arstar,armin,armax
18	format('angular radius of star =',f8.4,/,'angular radius of
     1  inner shell = ',f8.4,/,'angular radius of outer shell = ',
     2  f8.4,/, '(arcseconds)')
	write(1,5)
5	format(/)
c  Calculate total dust mass in envelope in solar masses
	xmd=(xml+xmll)*(rmax - rmin)/3.17d12/v0
	write(1,7) xmd
7	format('total dust mass = ',1pe10.3,'solar masses')
	write(1,5)
c  Calculate dust temperature at inner radius
	tg=tstar*(rstar/2./rmin)**0.4
	write(1,2) tg
2	format('grain temperature at inner radius = ',f7.1)
c  Generate frequency, wavelength and opacity grid
	do 300 j=1,np
	fh(j)=c*1.d4/xlam(j)
	freq(j)=fh(j)/1.d9
300	continue
c
c  Now modify grain opacity to be shallower at long wavelengths
c
	ind=30
	q0=q(ind)
	xlam0=xlam(ind)
	do 750 i=1,ind
	q(i)=q0*(xlam(i)/xlam0)**(-beta)
750	continue
c
c  Calculate optical depth versus frequency
c
	write(1,5)
	write(1,605)
605	format('calculate run of optical depth')
	write(1,5)
	const=xnd0*(1./rmin - 1./rmax)
	const=const*pi*ag**2
	constl=xnd0l*(1./rmin - 1./rmax)
	constl=constl*pi*agl**2
	do 600 i=1,np
	tau=const*q(i)+constl
	write(1,610) i,xlam(i),freq(i),tau,q(i)
600	continue
610	format(1x,i3,4(1pe10.3))
c Generate flux grid
	sum=0.
	sum1=0.
	do 140 i=1,np
	dum=dexp(h*fh(i)/xk/tstar) - 1.
	dum=c**2*dum
	dum=2.*pi*h*fh(i)**3/dum
	flux(i)=dum*4.*pi*rstar**2
	if(i.eq.1) then
	xint(i)=fh(i+1) - fh(i)
	go to 320
	end if
	if(i.eq.np) then
	xint(i)=fh(i)-fh(i-1)
	go to 320
	end if
	xint(i)=0.5*(fh(i+1)-fh(i-1))
320	continue
	sum=sum+flux(i)*xint(i)
	fjy(i)=flux(i)*flxcon
	sum1=sum1+fjy(i)*xint(i)
	fnudnu=fh(i)*fjy(i)
	write(4,8) i,freq(i),fjy(i),xlam(i),fnudnu
140	continue
	check=sum/solarl
	write(1,3) check
3	format('total luminosity = ',1pd10.3)
	sum1=sum1*4.*pi*(parsec*dist)**2/solarl
	sum1=sum1/xjansky
	write(1,3) sum1
	write(1,5)
	sum=0.
	dilute=4.*pi*rmin**2
	do 222 i=1,np
	flux(i)=flux(i)/dilute
222	continue
c  Calculate total dust mass
	xmd=0.
	tg=tg+20.
	xsmm=0.
	r1=rmin
	if(chi.lt.1.0) then
	chix=1.0
	else
	chix=chi
	end if
	adj=0.02/chix
	xindex=2.+adj
200	xindex=xindex-adj
	if(r2.eq.0.) go to 770
	fudge=4.*pi*r2**2*flxcon
	do 100 j=1,np
100	fjy(j)=fudge*flux(j)
770	continue
	if(xindex.lt.1.) xindex=1.
	dr=(adj)**xindex*r1
	r2=r1+dr
	if(r2.gt.rmax) r2=rmax
	if(r1.eq.r2) go to 2000
	ff=f0(4)/1.d9
	call sindex(fjy(14),fjy(15),freq(14),freq(15),
     &  ff,ss,gamma)
	fluxtrue(4)=ss
	fluxtrue(1)=fjy(77)
	fluxtrue(2)=fjy(51)
	fluxtrue(3)=fjy(28)
	do 220 j=1,4
	cor(j)=0.
220	fluxobs(j)=0.
	do 230 j=1,np
	do 235 k=1,4
	fluxobs(k)=fluxobs(k)+fjy(j)*filter(k,j)*xint(j)
c	cor(k)=cor(k)+fluxtrue(k)*f0(k)*filter(k,j)*xint(j)/fh(j)
 	cor(k)=cor(k)+f0(k)*filter(k,j)*xint(j)/fh(j)
c	cor(k)=cor(k)+filter(k,j)*xint(j)
235	continue
230	continue
	do 236 k=1,4
236	fluxobs(k)=fluxobs(k)/cor(k)
	write(12,12) r2,(fluxobs(k),k=1,4),(fluxtrue(k),k=1,4)
12	format(9(1pd10.3))
	ttt=r2/v0/1.d5/3.16d7
	col1=fluxobs(2)/fluxobs(1)
	col2=fluxobs(3)/fluxobs(2)
	col3=fluxobs(4)/fluxobs(3)
	col1=dlog10(col1)
	col2=dlog10(col2)
	col3=dlog10(col3)
	col1=2.5*col1
	col2=2.5*col2
	col3=2.5*col3
	s1=fluxobs(1)
	s2=fluxobs(2)
	s3=fluxobs(3)
	s4=fluxobs(4)
	f1=f0(1)
	f2=f0(2)
	f3=f0(3)
	f4=f0(4)
	f=1.d13
	call sindex(s1,s2,f1,f2,f,s,gamma1)
	call sindex(s2,s3,f2,f3,f,s,gamma2)
	call sindex(s3,s4,f3,f4,f,s,gamma3)
	write(13,13) r2,ttt,col1,col2,col3,gamma1,gamma2,gamma3
13	format(2(1pd10.3),1x,3(0pf8.3),3f6.2)
	xmd=xmd+4.*pi*(r2-r1)*(xnd0*xmg+xnd0l*xmgl)
	sum=0.
	suml=0.
c  Sum absorbed radiation per grain
	do 170 j=1,np
	sum=sum+flux(j)*q(j)*xint(j)
	suml=suml+flux(j)*xint(j)
170	continue
	xinput=sum
	xinputl=suml
2002	continue
c  Sum emitted radiation per grain
	quant0=8.*pi*h/c**2
	sum=0.
	do 400 j=1,np
	dum=h*fh(j)/xk/tg
	dum=(dexp(dum)-1.)
	dum=quant0*q(j)*fh(j)**3/dum
	sum=sum+dum*xint(j)
400	continue
	xoutput=sum
	diff=(xinput-xoutput)/xinput
	adiff=dabs(diff)
	if(adiff.lt.1.d-4) go to 2001
	if(diff.le.-1.d0) diff=-0.9
	tg=tg*(1.+diff)**0.2
	go to 2002
2001	continue
	tgl=(xinputl/4./sigma)**0.25
	write(1,6) r1,tg,tgl
6	format('shell distance from star = ',1pd10.3,2x,'grain 
     1  temperatures ',2(1x,0pf6.1))
	dens=xnd0/r1**2
	quant0=pi*ag**2*xnd0*(1./r1 - 1./r2)
	quant0l=pi*agl**2*xnd0l*(1./r1 - 1./r2)
	densl=xnd0l/r1**2
c
c  Now add in interstellar radiation field
c
c	go to 981
 	if(tg.lt.100.) then
 	dum=tg**4 + 18.**4
 	tg=dum**0.25
 	duml=tgl**4 + 18.**4
 	tgl=duml**0.25
 	write(1,6) r1,tg,tgl
 	write(1,5)
 	end if
  981   continue
	sum=0.
	rad=0.
	do 180 j=1,np
	quant=quant0*q(j)
	flux(j)=flux(j)*(1.-quant-quant0l)
	dum=h*fh(j)/xk/tg
	dum=(dexp(dum)-1.)*c**2
	dum=2.*pi*h*fh(j)**3/dum
 	dum=quant*dum*4.
	duml=h*fh(j)/xk/tgl
	duml=(dexp(duml)-1.)*c**2
	duml=2.*pi*h*fh(j)**3/duml
 	duml=quant0l*duml*4.
 	flux(j)=flux(j)+dum + duml
	flux(j)=flux(j)*(r1/r2)**2
	sum=sum+flux(j)*xint(j)
	rad=rad+xint(j)*dum
180	continue
	check=4.*pi*r2**2*sum/solarl
	write(1,3) check
	write(1,5)
	if(tg.gt.50.) then
	tg=tg*1.05
	else
	tg=tg
	end if
	r1=r2
	go to 200
2000	continue
	fudge=4.*pi*rmax**2*flxcon
	do 160 j=1,np
	flux(j)=fudge*flux(j)
	fnudnu=flux(j)*fh(j)
	write(3,8) j,freq(j),flux(j),xlam(j),fnudnu
8	format(1x,i3,4(1x,1pe10.3))
160	continue
	write(1,5)
	xmd=xmd/solarm
	write(1,7) xmd
1000	continue
	s1=fluxobs(1)
	s2=fluxobs(2)
	s3=fluxobs(3)
	s4=fluxobs(4)
	f1=f0(1)
	f2=f0(2)
	f3=f0(3)
	f4=f0(4)
	f=1.d13
	call sindex(s1,s2,f1,f2,f,s,gamma1)
	call sindex(s2,s3,f2,f3,f,s,gamma2)
	call sindex(s3,s4,f3,f4,f,s,gamma3)
	itstar=tstar
	idist=dist
	chi=dlog10(chi)
	do 3000 i=1,4
	f0(i)=f0(i)/1.d9
	write(8,3001) f0(i),fluxobs(i),fluxtrue(i)
3001	format(1x,3(2x,1pd10.4))
3000	continue
	stop
	end
c
	subroutine sindex(f1,f2,x1,x2,f,s,spind)
	implicit double precision (a-h, o-z)
	spind=(dlog10(f1)-dlog10(f2))/(dlog10(x1)-dlog10(x2))
	fudge=(f/x1)**spind
	s=f1*fudge
	return
	end