!*********************************************************************************
! Fortran code to calculate different species condensation flux boundaries.
!
! Ravi kumar Kopparapu, February 23, 2018
!
! Reference paper: Kopparapu et al. (2018), ApJ, arxiv: XXXXXXX
!*********************************************************************************
implicit none

real *8 flux(4),rp(6)
real *8 a(4), b(4), c(4), d(4), e(4), f(4)

integer i,j

! i = 1 ---> ZnS
! i = 2 ---> H2O
! i = 3 ---> CO2
! i = 4 ---> CH4

open(9,file='boundaries.dat')
open(10,file='boundaries_coeffs.dat')

a(1) = 0.0010338041*100.0d0
a(2) = 0.0017802416
a(3) = 0.0002947546
a(4) = -0.0000033096

b(1) = -0.0255230451 * 100.0d0
b(2) = -0.0443704003
b(3) = -0.0070583509
b(4) = 0.0000889715

c(1) = 0.1858822989 * 100.0d0
c(2) = 0.3302966853
c(3) = 0.0483928147
c(4) = -0.0007644190

d(1) = -0.4990171468 * 100.0d0
d(2) = -0.9299682851
d(3) = -0.1198359100
d(4) = 0.0026719183

e(1) = 0.5690844110 * 100.0d0
e(2) = 1.1366785108
e(3) = 0.1477297602
e(4) = -0.0038305535

f(1) = 1.6385396777 * 100.0d0
f(2) = 0.6255832476
f(3) = 0.2304769313
f(4) = 0.0048373922

rp(1) = 0.5
rp(2) = 1.0
rp(3) = 1.75
rp(4) = 3.5
rp(5) = 6.0
rp(6) = 14.3

100 Format('#',4x,'Rp(RE)',10x,'ZnS',10x,'H2O',12x,'CO2',10x,'CH4')
write(9,100)
do j = 1,6
  do i = 1,4
     flux(i) = (a(i)*rp(j)**5 + b(i)*rp(j)**4 + c(i)*rp(j)**3 + d(i)*rp(j)**2 +  e(i)*rp(j) + f(i))
  enddo
  write(9,200)rp(j),(flux(i),i=1,4)
enddo

write(10,*)'# The columns are a,b,c,d,e,f coefficients'
write(10,*)'# The rows are for different Rp/RE: 0.5, 1.0, 1.75, 3.5, 6, 14.3'
write(10,*)'# Note: The coefficients are applicable ONLY for the above radii'
write(10,*)'# Polynomial expression is:'
write(10,*)'#'
write(10,*)'#     a*Rp^5 + b*Rp^4 + c*Rp^3 + d*Rp^2 + e*Rp + f'
write(10,*)
write(10,300)(a(i),i=1,4)
write(10,300)(b(i),i=1,4)
write(10,300)(c(i),i=1,4)
write(10,300)(d(i),i=1,4)
write(10,300)(e(i),i=1,4)
write(10,300)(f(i),i=1,4)

200 format(1p8e14.5)
300 format(1p8e18.9)
close(9)
close(10)
end
