Difference between revisions of "Diels-Alder Reaction Experimental Design (VPLAN)"

c     RHS of the differential equations
 
      subroutine ffcn( t, x, f, p, q, rwh, iwh, iflag )
        implicit none
 
        real*8 x(*), f(*), p(*), q(*), rwh(*), t
        integer*4 iwh(*), iflag
 
        real*8 n1, n2, n3, n4
        real*8 na1, na2, na4
        real*8 fg, Temp, E , Rg , T1, Tc
        real*8 r1, mR
        real*8 kr1, kkat, Ckat, Ekat
        real*8 k1, lambda
        real*8 M1, M2, M3, M4
        real*8 dm
 
c       State variables
 
        n1 = x(1)
        n2 = x(2)
        n3 = x(3)
        n4 = x(4)
 
c       Control variables
 
        na1 = q(1)        
        na2 = q(2)
        na4 = q(3)
        Ckat = q(4)
 
c       Control function
 
c       DISCRETIZE1( Tc, rwh, iwh )
 
c       Parameters
 
        kr1 = p(1) * 1.0d-2 
        E = p(2) * 60000.0d+0
        k1 = p(3) * 0.10d+0 
        Ekat = p(4) * 40000.0d0
        lambda = p(5) * 0.25d+0
 
c       Molar masses (in kg/mol)
 
        M1 = 0.1362d+0
        M2 = 0.09806d+0
        M3 = M1 + M2
        M4 = 0.236d+0
 
        Temp = Tc + 273.0d+0
        Rg = 8.314d+0
        T1 = 293.0d+0
 
c       Reaction rates
 
        mR = n1*M1 + n2*M2 +n3*M3 + n4*M4
 
        kkat = kr1 * dexp( -E/Rg  * ( 1.0d+0/Temp - 1.0d+0/T1 ) )   
     &       + k1  * dexp( -Ekat/Rg *( 1.0d+0/Temp - 1.0d+0/T1 ) )
     &       * Ckat * dexp( -lambda * t  )
 
        r1 = kkat * n1 * n2 / mR
 
        f(1) = -r1           
        f(2) = -r1 
        f(3) = r1 
        f(4) = 0.0d0
 
      end

c     Dummyfunction for RHS of algebraic equations
 
      subroutine gfcn( t, x, g, p, q, rwh, iwh, iflag )
        implicit none
 
        real*8 x(*), g(*), p(*), q(*), rwh(*), t
        integer*4 iwh(*), iflag
 
        iflag=0
 
      end

c     Messfunktion
 
      subroutine mess3( t, x, h, p, q, rwh, iwh, iflag )
        implicit none
 
        real*8 t, x(*), h, p(*), q(*), rwh(*)
        integer*4 iwh(*), iflag
 
        real*8 M1, M2, M3, M4, mR
 
c       Berechnung der Reaktormasse
 
        M1 = 0.1362d+0
        M2 = 0.09806d+0
        M3 = M1 + M2
        M4 = 0.236d+0
 
        mR = M1*x(1) + M2*x(2) + M3*x(3) + M4*x(4)
 
c       Messwert: Massenprozent
 
        h = M3*x(3) * 100.0d+0/mR        
 
        iflag = 0
 
      end

c     Standardabweichung der Messfunktion
 
      subroutine sigma3( t, x, s, p, q, rwh, iwh, iflag )
        implicit none
 
        real*8 rwh(*)
        integer*4 iwh(*), iflag
 
        real*8 t, x(*), p(*), q(*)
        real*8 s
        real*8 h
 
        s = 1.0d+0
 
        iflag = 0
 
      end