MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
filmprop.f90
Go to the documentation of this file.
1 
7 module filmprop
8  use constants
9  use solidprop
10  implicit none
11  private
12  public scattwall,trextwall,absorwall,scattcoeffintwall,&
13  trextcoeffintwall,abscoeffintwall,filmconst
14 contains
15 !********************************BEGINNING*************************************
17 subroutine filmconst(lambda,theta,h,Rwin,Twin,Awin)
18  real(dp), intent(in) :: lambda
19  real(dp), intent(in) :: theta
20  real(dp), intent(in) :: h
21  real(dp), intent(out) :: Rwin
22  real(dp), intent(out) :: Twin
23  real(dp), intent(out) :: Awin
24  integer :: &
25  method=1 !1 is recommended, 3 and 4 should work; 2 is for no interference;
26  ! 1,3,4 give similar results
27  real(dp) :: &
28  theta2,& !refracted agle
29  rho12,rho23,& !reflectivities
30  tau12,tau23,& !transmissivities
31  r12,r23,& !reflection coefficients
32  t12,t23,& !transmission coefficients
33  dzeta,& !interference parameter
34  p,q
35  complex(dp) :: &
36  cn,beta,rc,tc,rc2,tc2
37  call optconst(lambda,n2,k2)
38  select case(method)
39  case(1)
40  !according to Modest
41  kappa2=4*pi*k2/lambda
42  ! theta2=asin(sin(theta)/n2)
43  ! rho12=(((cos(theta2)-n2*cos(theta))/(cos(theta2)+n2*cos(theta)))**2+&
44  ! ((cos(theta)-n2*cos(theta2))/(cos(theta)+n2*cos(theta2)))**2)/2
45  ! rho23=(((n2*cos(theta)-cos(theta2))/(n2*cos(theta)+cos(theta2)))**2+&
46  ! ((n2*cos(theta2)-cos(theta))/(n2*cos(theta2)+cos(theta)))**2)/2
47  ! write(*,*) rho23
48  ! stop
49  p=sqrt((sqrt((n2**2-k2**2-n1**2*sin(theta)**2)**2+4*n2**2*k2**2)+&
50  (n2**2-k2**2-n1**2*sin(theta)**2))/2)
51  q=sqrt((sqrt((n2**2-k2**2-n1**2*sin(theta)**2)**2+4*n2**2*k2**2)-&
52  (n2**2-k2**2-n1**2*sin(theta)**2))/2)
53  theta2=atan(sin(theta)*n1/p)
54  rho12=((n1*cos(theta)-p)**2+q**2)/((n1*cos(theta)+p)**2+q**2)
55  rho12=rho12*(1+((p-n1*sin(theta)*tan(theta))**2+q**2)/&
56  ((p+n1*sin(theta)*tan(theta))**2+q**2))/2
57  rho23=rho12
58  ! write(*,*) rho12
59  ! stop
60  r12=-sqrt(rho12)
61  r23=sqrt(rho23)
62  tau12=1-rho12
63  tau23=1-rho23
64  dzeta=4*pi*n2*h/lambda
65  rwin=(r12**2+2*r12*r23*exp(-kappa2*h)*cos(dzeta)+r23**2*exp(-2*kappa2*h))/&
66  (1+2*r12*r23*exp(-kappa2*h)*cos(dzeta)+r12**2*r23**2*exp(-2*kappa2*h))
67  twin=(tau12*tau23*exp(-kappa2*h))/&
68  (1+2*r12*r23*exp(-kappa2*h)*cos(dzeta)+r12**2*r23**2*exp(-2*kappa2*h))
69  awin=1-rwin-twin
70  case(2)
71  !no interference
72  kappa2=4*pi*k2/lambda
73  theta2=asin(sin(theta)/n2)
74  rho12=(((cos(theta2)-n2*cos(theta))/(cos(theta2)+n2*cos(theta)))**2+&
75  ((cos(theta)-n2*cos(theta2))/(cos(theta)+n2*cos(theta2)))**2)/2
76  rho23=(((n2*cos(theta)-cos(theta2))/(n2*cos(theta)+cos(theta2)))**2+&
77  ((n2*cos(theta2)-cos(theta))/(n2*cos(theta2)+cos(theta)))**2)/2
78  tau12=exp(-kappa2*h)
79  rwin=(rho12+(1-2*rho12)*rho23*tau12**2)/(1-rho12*rho23*tau12**2)
80  twin=(1-rho12)*(1-rho23)*tau12/(1-rho12*rho23*tau12**2)
81  awin=1-rwin-twin
82  case(3)
83  !Coquard (according to Ochsner), adapted
84  cn=cmplx(n2,-k2,kind=dp)
85  theta2=asin(sin(theta)/n2)
86  rho12=(((cos(theta2)-n2*cos(theta))/(cos(theta2)+n2*cos(theta)))**2+&
87  ((cos(theta)-n2*cos(theta2))/(cos(theta)+n2*cos(theta2)))**2)/2
88  rho23=(((n2*cos(theta)-cos(theta2))/(n2*cos(theta)+cos(theta2)))**2+&
89  ((n2*cos(theta2)-cos(theta))/(n2*cos(theta2)+cos(theta)))**2)/2
90  r12=-sqrt(rho12)
91  r23=sqrt(rho23)
92  tau12=1-r12**2
93  tau23=1-r23**2
94  t12=sqrt(tau12/n2)
95  t23=sqrt(tau23*n2)
96  beta=2*pi*cn*h*cos(theta2)/lambda
97  rc=(r12+r23*exp(-iu*2*beta))/(1+r12*r23*exp(-iu*2*beta))
98  tc=(t12*t23*exp(-iu*beta))/(1+r12*r23*exp(-iu*2*beta))
99  rwin=abs(rc)**2
100  twin=abs(tc)**2
101  awin=1-rwin-twin
102  case(4)
103  !Coquard (according to Dombrovsky)
104  cn=cmplx(n2,-k2,kind=dp)
105  theta2=asin(sin(theta)/n2)
106  r12=(cos(theta2)-n2*cos(theta))/(cos(theta2)+n2*cos(theta))
107  ! r12=(n2*cos(theta)-cos(theta2))/(n2*cos(theta)+cos(theta2)) !this gives
108  ! the same results
109  r23=-r12
110  t12=1-r12**2 !actually this is t12*t23
111  beta=2*pi*cn*h*cos(theta2)/lambda
112  rc=(r12+r23*exp(-iu*2*beta))/(1+r12*r23*exp(-iu*2*beta))
113  tc=(t12*exp(-iu*beta))/(1+r12*r23*exp(-iu*2*beta))
114  r12=(cos(theta)-n2*cos(theta2))/(cos(theta)+n2*cos(theta2))
115  ! r12=(n2*cos(theta2)-cos(theta))/(n2*cos(theta2)+cos(theta))
116  r23=-r12
117  t12=1-r12**2 !actually this is t12*t23
118  rc2=(r12+r23*exp(-iu*2*beta))/(1+r12*r23*exp(-iu*2*beta))
119  tc2=(t12*exp(-iu*beta))/(1+r12*r23*exp(-iu*2*beta))
120  rwin=(abs(rc)**2+abs(rc2)**2)/2
121  twin=(abs(tc)**2+abs(tc2)**2)/2
122  awin=1-rwin-twin
123  case default
124  stop 'unknown method for calculation of slab reflectivity'
125  end select
126 end subroutine filmconst
127 !***********************************END****************************************
128 
129 
130 !********************************BEGINNING*************************************
132 real(dp) function scattwall ( dw )
133  use quadpack
134  real(dp), intent(in) :: dw
135  real(dp) :: wtweight,dwold
136  !qags variables
137  real(dp) :: a !start point of integration
138  real(dp) :: abserr
139  real(dp) :: b !end point of integration
140  real(dp), parameter :: epsabs = 0.0e0_dp
141  real(dp), parameter :: epsrel = 1e-3_dp
142  integer :: ier
143  integer :: neval
144  real(dp) :: res
145  !end of qags variables
146  a=0
147  b=pi/2
148  wtweight=1/(sqrt(2*pi)*wsdev*dw)*exp(-(log(dw)-log(dwall))**2/(2*wsdev**2))
149  dwold=dwall
150  dwall=dw
151  call qags(scattcoeffintwall, a, b, epsabs, epsrel, res, abserr, neval, ier)
152  if (ier /= 0) then
153  write(*,*) 'qags returned',ier
154  write(*,*) 'wall thickness',dw
155  write(*,*) 'scattering coefficient of walls not calculated'
156  stop
157  endif
158  scattwall=res*(1-por)/dwall*wtweight
159  dwall=dwold
160 end function scattwall
161 !***********************************END****************************************
162 
163 
164 !********************************BEGINNING*************************************
166 real(dp) function trextwall ( dw )
167  use quadpack
168  real(dp), intent(in) :: dw
169  real(dp) :: wtweight,dwold
170  !qags variables
171  real(dp) :: a !start point of integration
172  real(dp) :: abserr
173  real(dp) :: b !end point of integration
174  real(dp), parameter :: epsabs = 0.0e0_dp
175  real(dp), parameter :: epsrel = 0.001e0_dp
176  integer :: ier
177  integer :: neval
178  real(dp) :: res
179  !end of qags variables
180  a=0
181  b=pi/2
182  wtweight=1/(sqrt(2*pi)*wsdev*dw)*exp(-(log(dw)-log(dwall))**2/(2*wsdev**2))
183  dwold=dwall
184  dwall=dw
185  call qags(trextcoeffintwall, a, b, epsabs, epsrel, res, abserr, neval, ier)
186  if (ier /= 0) then
187  write(*,*) 'qags returned',ier
188  write(*,*) 'wall thickness',dw
189  write(*,*) 'transport extinction coefficient of walls not calculated'
190  stop
191  endif
192  trextwall=res*(1-por)/dwall*wtweight
193  dwall=dwold
194 end function trextwall
195 !***********************************END****************************************
196 
197 
198 !********************************BEGINNING*************************************
200 real(dp) function absorwall ( dw )
201  use quadpack
202  real(dp), intent(in) :: dw
203  real(dp) :: wtweight,dwold
204  !qags variables
205  real(dp) :: a !start point of integration
206  real(dp) :: abserr
207  real(dp) :: b !end point of integration
208  real(dp), parameter :: epsabs = 0.0e0_dp
209  real(dp), parameter :: epsrel = 0.001e0_dp
210  integer :: ier
211  integer :: neval
212  real(dp) :: res
213  !end of qags variables
214  a=0
215  b=pi/2
216  wtweight=1/(sqrt(2*pi)*wsdev*dw)*exp(-(log(dw)-log(dwall))**2/(2*wsdev**2))
217  dwold=dwall
218  dwall=dw
219  call qags ( abscoeffintwall, a, b, epsabs, epsrel, res, abserr, neval, ier )
220  if (ier /= 0) then
221  write(*,*) 'qags returned',ier
222  write(*,*) 'wall thickness',dw
223  write(*,*) 'absorption coefficient of walls not calculated'
224  stop
225  endif
226  absorwall=res*(1-por)/dwall*wtweight
227  dwall=dwold
228 end function absorwall
229 !***********************************END****************************************
230 
231 
232 !********************************BEGINNING*************************************
234 real(dp) function scattcoeffintwall ( theta )
235  real(dp), intent(in) :: theta
236  real(dp) :: Rwin,Twin,Awin
237  call filmconst(lambda,theta,dwall,rwin,twin,awin)
238  scattcoeffintwall=rwin*sin(theta)*cos(theta)
239 end function scattcoeffintwall
240 !***********************************END****************************************
241 
242 
243 !********************************BEGINNING*************************************
245 real(dp) function trextcoeffintwall ( theta )
246  real(dp), intent(in) :: theta
247  real(dp) :: Rwin,Twin,Awin
248  call filmconst(lambda,theta,dwall,rwin,twin,awin)
249  trextcoeffintwall=(1-twin+rwin*cos(2*theta))*sin(theta)*cos(theta)
250 end function trextcoeffintwall
251 !***********************************END****************************************
252 
253 
254 !********************************BEGINNING*************************************
256 real(dp) function abscoeffintwall ( theta )
257  real(dp), intent(in) :: theta
258  real(dp) :: Rwin,Twin,Awin
259  call filmconst(lambda,theta,dwall,rwin,twin,awin)
260  abscoeffintwall=awin*sin(theta)*cos(theta)
261 end function abscoeffintwall
262 !***********************************END****************************************
263 end module filmprop
lambda
double lambda
Latent heat of blowing agent, J/kg.
Definition: experimentalInputs.h:118