📄 conf.f
字号:
mu4_variance=mu4_variance/DFLOAT((NSTEP-1)*Nsets)
conc_variance=conc_variance/DFLOAT((NSTEP-1)*Nsets)
DO 1370 I = 1, NTHETA
DO 1375 J = 1, Nu*(NSTEP-1)*Nsets, Nu
FTV(I,J)=F(J,I)/mu0_variance
FTV(I,J+1)=F(J+1,I)/mu1_variance
FTV(I,J+2)=F(J+2,I)/mu2_variance
FTV(I,J+3)=F(J+3,I)/mu3_variance
FTV(I,J+4)=F(J+4,I)/mu4_variance
FTV(I,J+5)=F(J+5,I)/conc_variance
1375 CONTINUE
1370 CONTINUE
DO 1400 I = 1, NTHETA
DO 1500 J = 1, NTHETA
FTVF(I,J) = 0.0D0
DO 1600 K = 1, Nu*(NSTEP-1)*Nsets
FTVF(I,J)=FTVF(I,J)+FTV(I,K)*F(K,J)
1600 CONTINUE
1500 CONTINUE
1400 CONTINUE
*Calculate determinant of FTVF,
*which is symmetric with dimensions NTHETA by NTHETA
detFTVF = (FTVF(1,3)*FTVF(2,4)-FTVF(1,4)*FTVF(2,3))**2.D0+
& (FTVF(1,3)*FTVF(3,4)-FTVF(3,3)*FTVF(1,4))*
& (FTVF(1,4)*FTVF(2,2)-FTVF(1,2)*FTVF(2,4))+
& (FTVF(2,3)*FTVF(3,4)-FTVF(3,3)*FTVF(2,4))*
& (FTVF(1,1)*FTVF(2,4)-FTVF(1,4)*FTVF(1,2))+
& (FTVF(1,3)*FTVF(4,4)-FTVF(3,4)*FTVF(1,4))*
& (FTVF(1,2)*FTVF(2,3)-FTVF(1,3)*FTVF(2,2))+
& (FTVF(2,3)*FTVF(4,4)-FTVF(3,4)*FTVF(2,4))*
& (FTVF(1,2)*FTVF(1,3)-FTVF(1,1)*FTVF(2,3))+
& (FTVF(3,3)*FTVF(4,4)-FTVF(3,4)*FTVF(3,4))*
& (FTVF(1,1)*FTVF(2,2)-FTVF(1,2)*FTVF(1,2))
PRINT*, "Determinant of FTVF = ", detFTVF
PRINT*, '-natural log (determinant of FTVF)= ',
& -DLOG(detFTVF)
PRINT*,'FTVF = '
WRITE(*,11)FTVF(1,1), FTVF(1,2),FTVF(1,3),FTVF(1,4)
WRITE(*,11)FTVF(2,1), FTVF(2,2),FTVF(2,3),FTVF(2,4)
WRITE(*,11)FTVF(3,1), FTVF(3,2),FTVF(3,3),FTVF(3,4)
WRITE(*,11)FTVF(4,1), FTVF(4,2),FTVF(4,3),FTVF(4,4)
pause
11 FORMAT(4(E24.16,2x))
CALL DLINRG(NTHETA,FTVF,NTHETA, inv_FTVF,NTHETA)
* From chi squared distribution
chi_squared=9.49D0
g_interval=DSQRT(chi_squared*inv_FTVF(1,1))
lnkg_interval=DSQRT(chi_squared*inv_FTVF(2,2))
b_interval=DSQRT(chi_squared*inv_FTVF(3,3))
lnkb_interval=DSQRT(chi_squared*inv_FTVF(4,4))
PRINT*
PRINT*,'Using chi-squared = ',chi_squared
PRINT*
PRINT*,'g = ', g, ' +/- ', g_interval
PRINT*,'ln kg = ', DLOG(kg), ' +/- ', lnkg_interval
PRINT*,' b = ', b, ' +/- ', b_interval
PRINT*,'ln kb = ', DLOG(kb), ' +/- ', lnkb_interval
PRINT*
PRINT*, 'mu0_variance = ', mu0_variance
PRINT*, 'mu1_variance = ', mu1_variance
PRINT*, 'mu2_variance = ', mu2_variance
PRINT*, 'mu3_variance = ', mu3_variance
PRINT*, 'mu4_variance = ', mu4_variance
PRINT*, 'conc_variance = ', conc_variance
PRINT*, 'mu0_sigma = ', DSQRT(mu0_variance)
PRINT*, 'mu1_sigma = ', DSQRT(mu1_variance)
PRINT*, 'mu2_sigma = ', DSQRT(mu2_variance)
PRINT*, 'mu3_sigma = ', DSQRT(mu3_variance)
PRINT*, 'mu4_sigma = ', DSQRT(mu4_variance)
PRINT*, 'conc_sigma = ', DSQRT(conc_variance)
OPEN(UNIT=50, FILE='FTVF.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=60, FILE='variance.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
OPEN(UNIT=70, FILE='interval.dat', FORM='FORMATTED',
& ACCESS='SEQUENTIAL', STATUS='UNKNOWN')
WRITE(50,11)FTVF(1,1), FTVF(1,2),FTVF(1,3),FTVF(1,4)
WRITE(50,11)FTVF(2,1), FTVF(2,2),FTVF(2,3),FTVF(2,4)
WRITE(50,11)FTVF(3,1), FTVF(3,2),FTVF(3,3),FTVF(3,4)
WRITE(50,11)FTVF(4,1), FTVF(4,2),FTVF(4,3),FTVF(4,4)
WRITE(60,12) mu0_variance
WRITE(60,12) mu1_variance
WRITE(60,12) mu2_variance
WRITE(60,12) mu3_variance
WRITE(60,12) mu4_variance
WRITE(60,12) conc_variance
WRITE(70,12) g_interval
WRITE(70,12) lnkg_interval
WRITE(70,12) b_interval
WRITE(70,12) lnkb_interval
CLOSE(UNIT=50)
CLOSE(UNIT=60)
CLOSE(UNIT=70)
CLOSE(UNIT=90)
12 FORMAT(E24.16)
print*, "done"
STOP
END
****************************************************************
SUBROUTINE FCN(NEQ,T,Y,Yprime)
INTEGER NEQ, I, J, K
INTEGER NTHETA, NU
REAL*8 T, Y(NEQ), Yprime(NEQ)
REAL*8 r0, alpha, mu00, kg, g, kb, b
REAL*8 Csat, birth, growth, Temp
REAL*8 JACOBIAN(6,6), W(6,4)
REAL*8 WPRIME(6,4), F_THETA(6,4)
COMMON /EXP_DATA/r0, alpha, mu00
COMMON /GROWTH_DATA/kg, g
COMMON /BIRTH_DATA/kb, b
NTHETA=4
NU=6
*Moments
* d(mu0(t))/dt, # of particles/g solvent/minute
Yprime(1)=birth(Y(6),Csat(Temp(T)),Y(4))
* d(mu1(t))/dt, microns/g solvent/minute
Yprime(2)=growth(Y(6),Csat(Temp(T)))*Y(1)+
& birth(Y(6),Csat(Temp(T)),Y(4))*r0
* d(mu2(t))/dt, microns^2/g solvent/minute
Yprime(3)=2.0D0*growth(Y(6),Csat(Temp(T)))*Y(2)+
& birth(Y(6),Csat(Temp(T)),Y(4))*r0**2
* d(mu3(t))/dt, microns^3/g solvent/minute
Yprime(4)=3.0D0*growth(Y(6),Csat(Temp(T)))*Y(3)+
& birth(Y(6),Csat(Temp(T)),Y(4))*r0**3
* d(mu4(t))/dt, microns^4/g solvent/minute
Yprime(5)=4.0D0*growth(Y(6),Csat(Temp(T)))*Y(4)+
& birth(Y(6),Csat(Temp(T)),Y(4))*r0**4
*Mass balance
* d(conc(t))/dt, g solute/g solvent/minute
Yprime(6)=-alpha*(3.0D0*growth(Y(6),Csat(Temp(T)))*Y(3)+
& birth(Y(6),Csat(Temp(T)),Y(4))*r0**3)
*Moments for seed crystals only
* d(mu'1(t)/dt), microns/g solvent/minute
Yprime(7)=growth(Y(6),Csat(Temp(T)))*mu00
* d(mu'2(t)/dt), microns^2/g solvent/minute
Yprime(8)=2.0D0*growth(Y(6),Csat(Temp(T)))*Y(7)
* d(mu'2]3(t)/dt), microns^3/g solvent/minute
Yprime(9)=3.0D0*growth(Y(6),Csat(Temp(T)))*Y(8)
DO 203 I = 1 , NU
DO 204 J = 1 , NU
JACOBIAN(I,J)=0.0D0
204 CONTINUE
203 CONTINUE
*The following equations assume r0 = 0 (nucleation crystal size)
JACOBIAN(1,4)=birth(Y(6),Csat(Temp(T)),Y(4))/Y(4)
JACOBIAN(1,6)=b*birth(Y(6),Csat(Temp(T)),Y(4))/
& (Y(6)-Csat(Temp(T)))
JACOBIAN(2,1)=growth(Y(6),Csat(Temp(T)))
JACOBIAN(2,6)=g*growth(Y(6),Csat(Temp(T)))*Y(1)/
& (Y(6)-Csat(Temp(T)))
JACOBIAN(3,2)=2.0D0*growth(Y(6),Csat(Temp(T)))
JACOBIAN(3,6)=2.0D0*g*growth(Y(6),Csat(Temp(T)))*Y(2)/
& (Y(6)-Csat(Temp(T)))
JACOBIAN(4,3)=3*growth(Y(6),Csat(Temp(T)))
JACOBIAN(4,6)=3.0D0*g*growth(Y(6),Csat(Temp(T)))*Y(3)/
& (Y(6)-Csat(Temp(T)))
JACOBIAN(5,4)=4*growth(Y(6),Csat(Temp(T)))
JACOBIAN(5,6)=4.0D0*g*growth(Y(6),Csat(Temp(T)))*Y(4)/
& (Y(6)-Csat(Temp(T)))
JACOBIAN(6,3)=-3.0D0*alpha*growth(Y(6),Csat(Temp(T)))
JACOBIAN(6,6)=-3.0D0*alpha*g*growth(Y(6),Csat(Temp(T)))*Y(3)
& /(Y(6)-Csat(Temp(T)))
K=10
DO 213 I =1 , NU
DO 214 J = 1, NTHETA
W(I,J)=Y(K)
K=K+1
214 CONTINUE
213 CONTINUE
DO 223 I = 1 , NU
DO 224 J = 1 , NTHETA
F_THETA(I,J)=0.0D0
224 CONTINUE
223 CONTINUE
*Again, the following equations assume r0 = 0
F_THETA(1,3)=birth(Y(6),Csat(Temp(T)),Y(4))*
& DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
F_THETA(1,4)=birth(Y(6), Csat(Temp(T)),Y(4))
F_THETA(2,1)=Y(1)*growth(Y(6),Csat(Temp(T)))*
& DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
F_THETA(2,2)=Y(1)*growth(Y(6),Csat(Temp(T)))
F_THETA(3,1)=2.0D0*Y(2)*growth(Y(6),Csat(Temp(T)))*
& DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
F_THETA(3,2)=2.0D0*Y(2)*growth(Y(6),Csat(Temp(T)))
F_THETA(4,1)=3.0D0*Y(3)*growth(Y(6),Csat(Temp(T)))*
& DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
F_THETA(4,2)=3.0D0*Y(3)*growth(Y(6),Csat(Temp(T)))
F_THETA(5,1)=4.0D0*Y(4)*growth(Y(6),Csat(Temp(T)))*
& DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
F_THETA(5,2)=4.0D0*Y(4)*growth(Y(6),Csat(Temp(T)))
F_THETA(6,1)=-3.0D0*alpha*Y(3)*growth(Y(6),Csat(Temp(T)))*
& DLOG(DABS((Y(6)-Csat(Temp(T)))/Csat(Temp(T))))
F_THETA(6,2)=-3.0D0*alpha*Y(3)*growth(Y(6),Csat(Temp(T)))
*Sensitivity Equation: W' = J W + d f/d theta
*where
* d u/ d t = f(t,u;theta) (solved above)
* W = d u /d theta (unknown)
* J = d f / du (Jacobian matrix)
DO 101 I = 1,NU
DO 102 J = 1,NTHETA
WPRIME(I,J)=F_THETA(I,J)
102 CONTINUE
101 CONTINUE
DO 521 I = 1,NU
DO 522 J = 1,NTHETA
DO 523 K = 1,NU
WPRIME(I,J)=WPRIME(I,J)+JACOBIAN(I,K)*W(K,J)
523 CONTINUE
522 CONTINUE
521 CONTINUE
K=10
DO 243 I = 1 , NU
DO 244 J = 1 , NTHETA
Yprime(K)=WPRIME(I,J)
K=K+1
244 CONTINUE
243 CONTINUE
RETURN
END
****************************************************************
SUBROUTINE FCNJ(NEQ,T,Y,DYPDY)
INTEGER NEQ
REAL*8 T,Y(NEQ),DYPDY(NEQ,*)
RETURN
END
****************************************************************
REAL*8 FUNCTION Temp(time)
* Crystallizer temperature setpoint profile for the
* simulation of a batch cooling crystallizer
*
* input: time - minutes
* output: Temp - temperature in degrees Centigrade
REAL*8 time
* Temp=(28.0D0-32.0D0)/160.0D0*time+32.0D0
Temp=(28.0D0-32.0D0)*(1.0D0-DEXP(-time/310.0D0))+32.0D0
RETURN
END
*************************************************************
REAL*8 FUNCTION Csat(T)
* saturation concentration for the simulation of a
* cooling batch crystallizer (potassium nitrate-water)
* system, from Appendix C in Miller
*
* input: T - temperature (20-40 degree Centigrade)
* output: Csat - saturation concentration
* (g KNO3/g water)
REAL*8 T
Csat=0.1286D0+0.00588D0*T+0.0001721D0*T**2
RETURN
END
****************************************************************
REAL*8 FUNCTION heat(conc)
* heat capacity for the simulation of a cooling
* batch crystallizer (potassium nitrate-water
* system), from Appendix C in Miller
*
* input: conc - soluction concentration
* output: heat - heat capacity, calories/g solvent/ degree C
REAL*8 conc
heat=0.9971-1.06*conc/(1+conc)+1.007*(conc/(1+conc))**2
RETURN
END
*****************************************************************
REAL*8 FUNCTION enthalpy(conc)
* enthalpy of crystallization for the simulation of a
* cooling batch crystallizer (potassium nitrate-water
* system), from Appendix C in Miller
*
* input: conc - solute concentration
* output: enthalpy - enthalpy of crystallization
* (calories/g solvent)
REAL*8 conc
enthalpy=-85.75+92.82*conc-99.97*conc**2
RETURN
END
************************************************************
REAL*8 FUNCTION growth(conc, concs)
* growth rate for the simulation of a cooling batch
* crystallizer
*
* arguments: conc - solute concentration
* concs - saturation concentration
* non-argument input: kg, g kinetic rate parameters
* output: growth - growth rate
REAL*8 conc, concs, kg, g
COMMON /GROWTH_DATA/kg, g
if ((conc-concs) .gt. 0) then
growth=kg*((conc-concs)/concs)**g
elseif ((conc-concs) .eq. 0) then
growth = 0.0D0
else
growth = -1.0D0*kg*(-1.0D0*(conc-concs)/concs)**g
endif
RETURN
END
************************************************************
REAL*8 FUNCTION birth(conc, concs,m3)
* birthth rate for the simulation of a cooling batch
* crystallizer
*
* arguments: conc - solute concentration
* concs - saturation concentration
* m3 - 3rd moment
* non-argument input: kb, b kinetic rate parameters
* output: birth - birth rate
REAL*8 conc, concs, m3, kb, b
COMMON /BIRTH_DATA/kb, b
if ((conc-concs) .gt. 0) then
birth = kb*(((conc-concs)/concs)**b)*m3*(1.0D-4)**3.0D0
elseif ((conc-concs) .eq. 0) then
birth = 0.0D0
else
birth = -1.0D0*kb*((-1.0D0*(conc-concs)/concs)**b)
& *m3*(1.0D-4)**3.0D0
endif
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -