psins.r

r软件 另一款可以计算核估计的软件包 需安装r软件

    psi <- -3*(2*a13*a23^2 + a13*a22*a33 + 2*a12*a23*a33)
  if ((r[1]==1) & (r[2]==2) & (r[3]==2) & (r[4]==1))
    psi <- -2*a14*a23^2 - 4*a13*a23*a24 - a14*a22*a33 - 2*a12*a24*a33 -
      2*a13*a22*a34 -  4*a12*a23*a34
  if ((r[1]==1) & (r[2]==2) & (r[3]==1) & (r[4]==2))
    psi <- -4*a14*a23*a24 - 2*a13*a24^2 - 2*a14*a22*a34 - 4*a12*a24*a34 - a13*a22*a44 - 2*a12*a23*a44
  if ((r[1]==1) & (r[2]==2) & (r[3]==0) & (r[4]==3))
    psi <- -3*(2*a14*a24^2 + a14*a22*a44 + 2*a12*a24*a44)
  if ((r[1]==1) & (r[2]==1) & (r[3]==4) & (r[4]==0))
    psi <- -3*a33*(4*a13*a23 + a12*a33)
  if ((r[1]==1) & (r[2]==1) & (r[3]==3) & (r[4]==1))
    psi <- -3*(a14*a23*a33 + a13*a24*a33 + 2*a13*a23*a34 + a12*a33*a34)
  if ((r[1]==1) & (r[2]==1) & (r[3]==2) & (r[4]==2))
    psi <- -2*a14*a24*a33 - 4*a14*a23*a34 - 4*a13*a24*a34 - 2*a12*a34^2 - 
      2*a13*a23*a44 - a12*a33*a44
  if ((r[1]==1) & (r[2]==1) & (r[3]==1) & (r[4]==3))
    psi <- -3*(2*a14*a24*a34 + a14*a23*a44 + a13*a24*a44 + a12*a34*a44)
  if ((r[1]==1) & (r[2]==1) & (r[3]==0) & (r[4]==4))
    psi <- -3*a44*(4*a14*a24 + a12*a44)
  if ((r[1]==1) & (r[2]==0) & (r[3]==5) & (r[4]==0))
    psi <- -15*a13*a33^2
  if ((r[1]==1) & (r[2]==0) & (r[3]==4) & (r[4]==1))
    psi <- -3*a33*(a14*a33 + 4*a13*a34)
  if ((r[1]==1) & (r[2]==0) & (r[3]==3) & (r[4]==2))
    psi <- -3*(2*a14*a33*a34 + 2*a13*a34^2 + a13*a33*a44)
  if ((r[1]==1) & (r[2]==0) & (r[3]==2) & (r[4]==3))
    psi <- -3*(2*a14*a34^2 + a14*a33*a44 + 2*a13*a34*a44)
  if ((r[1]==1) & (r[2]==0) & (r[3]==1) & (r[4]==4))
    psi <- -3*a44*(4*a14*a34 + a13*a44)
  if ((r[1]==1) & (r[2]==0) & (r[3]==0) & (r[4]==5))
    psi <- -15*a14*a44^2
  if ((r[1]==0) & (r[2]==6) & (r[3]==0) & (r[4]==0))
    psi <- -15*a22^3
  if ((r[1]==0) & (r[2]==5) & (r[3]==1) & (r[4]==0))
    psi <- -15*a22^2*a23
  if ((r[1]==0) & (r[2]==5) & (r[3]==0) & (r[4]==1))
    psi <- -15*a22^2*a24
  if ((r[1]==0) & (r[2]==4) & (r[3]==2) & (r[4]==0))
    psi <- -3*a22*(4*a23^2 + a22*a33)
  if ((r[1]==0) & (r[2]==4) & (r[3]==1) & (r[4]==1))
    psi <- -3*a22*(4*a23*a24 + a22*a34)
  if ((r[1]==0) & (r[2]==4) & (r[3]==0) & (r[4]==2))
    psi <- -3*a22*(4*a24^2 + a22*a44)
  if ((r[1]==0) & (r[2]==3) & (r[3]==3) & (r[4]==0))
    psi <- -6*a23^3 - 9*a22*a23*a33
  if ((r[1]==0) & (r[2]==3) & (r[3]==2) & (r[4]==1))
    psi <- -3*(2*a23^2*a24 + a22*a24*a33 + 2*a22*a23*a34)
  if ((r[1]==0) & (r[2]==3) & (r[3]==1) & (r[4]==2))
    psi <- -3*(2*a23*a24^2 + 2*a22*a24*a34 + a22*a23*a44)
  if ((r[1]==0) & (r[2]==3) & (r[3]==0) & (r[4]==3))
   psi <- -6*a24^3 - 9*a22*a24*a44
  if ((r[1]==0) & (r[2]==2) & (r[3]==4) & (r[4]==0))
    psi <- -3*a33*(4*a23^2 + a22*a33)
  if ((r[1]==0) & (r[2]==2) & (r[3]==3) & (r[4]==1))
    psi <- -3*(2*a23*a24*a33 + 2*a23^2*a34 + a22*a33*a34)
  if ((r[1]==0) & (r[2]==2) & (r[3]==2) & (r[4]==2))
    psi <- -2*a24^2*a33 - 8*a23*a24*a34 - 2*a22*a34^2 - 2*a23^2*a44 - a22*a33*a44
  if ((r[1]==0) & (r[2]==2) & (r[3]==1) & (r[4]==3))
    psi <- -3*(2*a24^2*a34 + 2*a23*a24*a44 + a22*a34*a44)
  if ((r[1]==0) & (r[2]==2) & (r[3]==0) & (r[4]==4))
    psi <- -3*a44*(4*a24^2 + a22*a44)
  if ((r[1]==0) & (r[2]==1) & (r[3]==5) & (r[4]==0))
    psi <- -15*a23*a33^2
  if ((r[1]==0) & (r[2]==1) & (r[3]==4) & (r[4]==1))
    psi <- -3*a33*(a24*a33 + 4*a23*a34)
  if ((r[1]==0) & (r[2]==1) & (r[3]==3) & (r[4]==2))
    psi <- -3*(2*a24*a33*a34 + 2*a23*a34^2 + a23*a33*a44)
  if ((r[1]==0) & (r[2]==1) & (r[3]==2) & (r[4]==3))
    psi <- -3*(2*a24*a34^2 + a24*a33*a44 + 2*a23*a34*a44)
  if ((r[1]==0) & (r[2]==1) & (r[3]==1) & (r[4]==4))
    psi <- -3*a44*(4*a24*a34 + a23*a44)
  if ((r[1]==0) & (r[2]==1) & (r[3]==0) & (r[4]==5))
    psi <- -15*a24*a44^2
  if ((r[1]==0) & (r[2]==0) & (r[3]==6) & (r[4]==0))
    psi <- -15*a33^3
  if ((r[1]==0) & (r[2]==0) & (r[3]==5) & (r[4]==1))
    psi <- -15*a33^2*a34
  if ((r[1]==0) & (r[2]==0) & (r[3]==4) & (r[4]==2))
    psi <- -3*a33*(4*a34^2 + a33*a44)
  if ((r[1]==0) & (r[2]==0) & (r[3]==3) & (r[4]==3))
    psi <- -6*a34^3 - 9*a33*a34*a44
  if ((r[1]==0) & (r[2]==0) & (r[3]==2) & (r[4]==4))
    psi <- -3*a44*(4*a34^2 + a33*a44)
  if ((r[1]==0) & (r[2]==0) & (r[3]==1) & (r[4]==5))
    psi <- -15*a34*a44^2
  if ((r[1]==0) & (r[2]==0) & (r[3]==0) & (r[4]==6))
    psi <- -15*a44^3

   return(psi)
}

psins.4d.part8 <- function(r, a11, a22, a33, a44, a12, a13, a14, a23, a24, a34)
{

  if ((r[1]==8) & (r[2]==0) & (r[3]==0) & (r[4]==0))
    psi <-105*a11^4
  if ((r[1]==7) & (r[2]==1) & (r[3]==0) & (r[4]==0))
    psi <- 105*a11^3*a12
  if ((r[1]==7) & (r[2]==0) & (r[3]==1) & (r[4]==0))
    psi <-105*a11^3*a13
  if ((r[1]==7) & (r[2]==0) & (r[3]==0) & (r[4]==1))
    psi <- 105*a11^3*a14
  if ((r[1]==6) & (r[2]==2) & (r[3]==0) & (r[4]==0))
    psi <- 15*a11^2*(6*a12^2 + a11*a22)
  if ((r[1]==6) & (r[2]==1) & (r[3]==1) & (r[4]==0))
    psi <- 15*a11^2*(6*a12*a13 + a11*a23)
  if ((r[1]==6) & (r[2]==1) & (r[3]==0) & (r[4]==1))
    psi <- 15*a11^2*(6*a12*a14 + a11*a24)
  if ((r[1]==6) & (r[2]==0) & (r[3]==2) & (r[4]==0))
    psi <- 15*a11^2*(6*a13^2 + a11*a33)
  if ((r[1]==6) & (r[2]==0) & (r[3]==1) & (r[4]==1))
    psi <- 15*a11^2*(6*a13*a14 + a11*a34)
  if ((r[1]==6) & (r[2]==0) & (r[3]==0) & (r[4]==2))
    psi <- 15*a11^2*(6*a14^2 + a11*a44)
  if ((r[1]==5) & (r[2]==3) & (r[3]==0) & (r[4]==0))
    psi <- 15*a11*a12*(4*a12^2 + 3*a11*a22)
  if ((r[1]==5) & (r[2]==2) & (r[3]==1) & (r[4]==0))
    psi <-15*a11*(4*a12^2*a13 + a11*a13*a22 + 2*a11*a12*a23)
  if ((r[1]==5) & (r[2]==2) & (r[3]==0) & (r[4]==1))
    psi <-15*a11*(4*a12^2*a14 + a11*a14*a22 + 2*a11*a12*a24)
  if ((r[1]==5) & (r[2]==1) & (r[3]==2) & (r[4]==0))
    psi <-15*a11*(4*a12*a13^2 + 2*a11*a13*a23 + a11*a12*a33)
  if ((r[1]==5) & (r[2]==1) & (r[3]==1) & (r[4]==1))
    psi <-15*a11*(4*a12*a13*a14 + a11*a14*a23 + a11*a13*a24 + a11*a12*a34)
  if ((r[1]==5) & (r[2]==1) & (r[3]==0) & (r[4]==2))
    psi <-15*a11*(4*a12*a14^2 + 2*a11*a14*a24 + a11*a12*a44)
  if ((r[1]==5) & (r[2]==0) & (r[3]==3) & (r[4]==0))
    psi <-15*a11*a13*(4*a13^2 + 3*a11*a33)
  if ((r[1]==5) & (r[2]==0) & (r[3]==2) & (r[4]==1))
    psi <-15*a11*(4*a13^2*a14 + a11*a14*a33 + 2*a11*a13*a34)
  if ((r[1]==5) & (r[2]==0) & (r[3]==1) & (r[4]==2))
    psi <-15*a11*(4*a13*a14^2 + 2*a11*a14*a34 + a11*a13*a44)
  if ((r[1]==5) & (r[2]==0) & (r[3]==0) & (r[4]==3))
    psi <-15*a11*a14*(4*a14^2 + 3*a11*a44)
  if ((r[1]==4) & (r[2]==4) & (r[3]==0) & (r[4]==0))
    psi <-24*a12^4 + 72*a11*a12^2*a22 + 9*a11^2*a22^2
  if ((r[1]==4) & (r[2]==3) & (r[3]==1) & (r[4]==0))
    psi <-24*a12^3*a13 + 36*a11*a12*a13*a22 + 36*a11*a12^2*a23 + 9*a11^2*a22*a23
  if ((r[1]==4) & (r[2]==3) & (r[3]==0) & (r[4]==1))
    psi <-24*a12^3*a14 + 36*a11*a12*a14*a22 + 36*a11*a12^2*a24 + 9*a11^2*a22*a24
  if ((r[1]==4) & (r[2]==2) & (r[3]==2) & (r[4]==0))
    psi <-
      3*(16*a11*a12*a13*a23 + 4*a12^2*(2*a13^2 + a11*a33) + 
         a11*(4*a13^2*a22 + 2*a11*a23^2 + a11*a22*a33))
  if ((r[1]==4) & (r[2]==2) & (r[3]==1) & (r[4]==1))
    psi <- 3*(8*a11*a12*(a14*a23 + a13*a24) + 4*a12^2*(2*a13*a14 + a11*a34) + 
              a11*(4*a13*a14*a22 + 2*a11*a23*a24 + a11*a22*a34))
  if ((r[1]==4) & (r[2]==2) & (r[3]==0) & (r[4]==2))
    psi <-
      3*(16*a11*a12*a14*a24 + 4*a12^2*(2*a14^2 + a11*a44) + 
         a11*(4*a14^2*a22 + 2*a11*a24^2 + a11*a22*a44))
  if ((r[1]==4) & (r[2]==1) & (r[3]==3) & (r[4]==0))
    psi <-
      9*a11*a23*(4*a13^2 + a11*a33) + 12*a12*(2*a13^3 + 3*a11*a13*a33)
  if ((r[1]==4) & (r[2]==1) & (r[3]==2) & (r[4]==1))
    psi <-
      3*(4*a12*(2*a13^2*a14 + a11*a14*a33 + 2*a11*a13*a34) + 
  a11*(8*a13*a14*a23 + 4*a13^2*a24 + a11*a24*a33 + 2*a11*a23*a34))
  if ((r[1]==4) & (r[2]==1) & (r[3]==1) & (r[4]==2))
    psi <-
      3*(4*a12*(2*a13*a14^2 + 2*a11*a14*a34 + a11*a13*a44) + 
         a11*(4*a14^2*a23 + 8*a13*a14*a24 + 2*a11*a24*a34 + a11*a23*a44))
  if ((r[1]==4) & (r[2]==1) & (r[3]==0) & (r[4]==3))
    psi <-
      9*a11*a24*(4*a14^2 + a11*a44) + 12*a12*(2*a14^3 + 3*a11*a14*a44)
  if ((r[1]==4) & (r[2]==0) & (r[3]==4) & (r[4]==0))
    psi <-
      24*a13^4 + 72*a11*a13^2*a33 + 9*a11^2*a33^2
  if ((r[1]==4) & (r[2]==0) & (r[3]==3) & (r[4]==1))
  psi <-
    24*a13^3*a14 + 36*a11*a13*a14*a33 + 36*a11*a13^2*a34 + 9*a11^2*a33*a34
  if ((r[1]==4) & (r[2]==0) & (r[3]==2) & (r[4]==2))
    psi <-
      3*(16*a11*a13*a14*a34 + 4*a13^2*(2*a14^2 + a11*a44) + 
         a11*(4*a14^2*a33 + 2*a11*a34^2 + a11*a33*a44))
  if ((r[1]==4) & (r[2]==0) & (r[3]==1) & (r[4]==3))
    psi <-
      9*a11*a34*(4*a14^2 + a11*a44) + 12*a13*(2*a14^3 + 3*a11*a14*a44)
  if ((r[1]==4) & (r[2]==0) & (r[3]==0) & (r[4]==4))
    psi <-
      24*a14^4 + 72*a11*a14^2*a44 + 9*a11^2*a44^2
  if ((r[1]==3) & (r[2]==5) & (r[3]==0) & (r[4]==0))
    psi <-
      15*a12*a22*(4*a12^2 + 3*a11*a22)
  if ((r[1]==3) & (r[2]==4) & (r[3]==1) & (r[4]==0))
    psi <-
      36*a12^2*a13*a22 + 9*a11*a13*a22^2 + 24*a12^3*a23 + 36*a11*a12*a22*a23
  if ((r[1]==3) & (r[2]==4) & (r[3]==0) & (r[4]==1))
    psi <-
      36*a12^2*a14*a22 + 9*a11*a14*a22^2 + 24*a12^3*a24 + 36*a11*a12*a22*a24
  if ((r[1]==3) & (r[2]==3) & (r[3]==2) & (r[4]==0))
    psi <-
      3*(12*a12^2*a13*a23 + 6*a11*a13*a22*a23 + 2*a12^3*a33 + 
         3*a12*(2*a13^2*a22 + 2*a11*a23^2 + a11*a22*a33))
  if ((r[1]==3) & (r[2]==3) & (r[3]==1) & (r[4]==1))
    psi <-
      3*(6*a12^2*(a14*a23 + a13*a24) + 3*a11*a22*(a14*a23 + a13*a24) + 
         2*a12^3*a34 + 3*a12*(2*a13*a14*a22 + 2*a11*a23*a24 + a11*a22*a34))
  if ((r[1]==3) & (r[2]==3) & (r[3]==0) & (r[4]==2))
    psi <-
      3*(12*a12^2*a14*a24 + 6*a11*a14*a22*a24 + 2*a12^3*a44 + 
         3*a12*(2*a14^2*a22 + 2*a11*a24^2 + a11*a22*a44))
  if ((r[1]==3) & (r[2]==2) & (r[3]==3) & (r[4]==0))
    psi <-
      3*(2*a13^3*a22 + 12*a12*a13^2*a23 + 6*a11*a12*a23*a33 + 
         3*a13*(2*a11*a23^2 + 2*a12^2*a33 + a11*a22*a33))
  if ((r[1]==3) & (r[2]==2) & (r[3]==2) & (r[4]==1))
    psi <-
      3*(2*a13^2*(a14*a22 + 2*a12*a24) + 2*a12^2*a14*a33 + 
         2*a13*(4*a12*a14*a23 + 2*a11*a23*a24 + 2*a12^2*a34 + a11*a22*a34) + 
         a11*(2*a14*a23^2 + a14*a22*a33 + 2*a12*a24*a33 + 4*a12*a23*a34))
  if ((r[1]==3) & (r[2]==2) & (r[3]==1) & (r[4]==2))
    psi <-
      3*(a13*(2*a14^2*a22 + 8*a12*a14*a24 + 2*a11*a24^2 + 2*a12^2*a44 + 
              a11*a22*a44) + 2*(2*a12^2*a14*a34 + a11*a14*(2*a23*a24 + a22*a34) + 
                              a12*(2*a14^2*a23 + 2*a11*a24*a34 + a11*a23*a44)))
  if ((r[1]==3) & (r[2]==2) & (r[3]==0) & (r[4]==3))
    psi <-
      3*(2*a14^3*a22 + 12*a12*a14^2*a24 + 6*a11*a12*a24*a44 + 
         3*a14*(2*a11*a24^2 + 2*a12^2*a44 + a11*a22*a44))
  if ((r[1]==3) & (r[2]==1) & (r[3]==4) & (r[4]==0))
    psi <-
      24*a13^3*a23 + 36*a12*a13^2*a33 + 36*a11*a13*a23*a33 + 9*a11*a12*a33^2
  if ((r[1]==3) & (r[2]==1) & (r[3]==3) & (r[4]==1))
    psi <-
      3*(2*a13^3*a24 + 6*a13^2*(a14*a23 + a12*a34) + 
         3*a11*a33*(a14*a23 + a12*a34) + 3*a13*(2*a12*a14*a33 + a11*a24*a33 + 
                                                2*a11*a23*a34))
  if ((r[1]==3) & (r[2]==1) & (r[3]==2) & (r[4]==2))
    psi <-
      3*(2*a11*a14*(a24*a33 + 2*a23*a34) + 2*a13^2*(2*a14*a24 + a12*a44) + 
         2*a13*(2*a14^2*a23 + 4*a12*a14*a34 + 2*a11*a24*a34 + a11*a23*a44) + 
         a12*(2*a14^2*a33 + 2*a11*a34^2 + a11*a33*a44))
  if ((r[1]==3) & (r[2]==1) & (r[3]==1) & (r[4]==3))
    psi <-
      3*(2*a14^3*a23 + 6*a14^2*(a13*a24 + a12*a34) + 3*a11*(a13*a24 + a12*a34)*
         a44 + 3*a14*(2*a11*a24*a34 + 2*a12*a13*a44 + a11*a23*a44))
  if ((r[1]==3) & (r[2]==1) & (r[3]==0) & (r[4]==4))
    psi <-
      24*a14^3*a24 + 36*a12*a14^2*a44 + 36*a11*a14*a24*a44 + 9*a11*a12*a44^2
  if ((r[1]==3) & (r[2]==0) & (r[3]==5) & (r[4]==0))
    psi <-
      15*a13*a33*(4*a13^2 + 3*a11*a33)
  if ((r[1]==3) & (r[2]==0) & (r[3]==4) & (r[4]==1))
    psi <-
      36*a13^2*a14*a33 + 9*a11*a14*a33^2 + 24*a13^3*a34 + 36*a11*a13*a33*a34
  if ((r[1]==3) & (r[2]==0) & (r[3]==3) & (r[4]==2))
    psi <-
      3*(12*a13^2*a14*a34 + 6*a11*a14*a33*a34 + 2*a13^3*a44 + 
         3*a13*(2*a14^2*a33 + 2*a11*a34^2 + a11*a33*a44))
  if ((r[1]==3) & (r[2]==0) & (r[3]==2) & (r[4]==3))
    psi <-
      3*(2*a14^3*a33 + 12*a13*a14^2*a34 + 6*a11*a13*a34*a44 + 
         3*a14*(2*a11*a34^2 + 2*a13^2*a44 + a11*a33*a44))
  if ((r[1]==3) & (r[2]==0) & (r[3]==1) & (r[4]==4))
    psi <-
      24*a14^3*a34 + 36*a13*a14^2*a44 + 36*a11*a14*a34*a44 + 9*a11*a13*a44^2
  if ((r[1]==3) & (r[2]==0) & (r[3]==0) & (r[4]==5))
    psi <-
      15*a14*a44*(4*a14^2 + 3*a11*a44)
  if ((r[1]==2) & (r[2]==6) & (r[3]==0) & (r[4]==0))
    psi <-
      15*a22^2*(6*a12^2 + a11*a22)
  if ((r[1]==2) & (r[2]==5) & (r[3]==1) & (r[4]==0))
    psi <-
      15*a22*(2*a12*a13*a22 + 4*a12^2*a23 + a11*a22*a23)
  if ((r[1]==2) & (r[2]==5) & (r[3]==0) & (r[4]==1))
    psi <-
      15*a22*(2*a12*a14*a22 + 4*a12^2*a24 + a11*a22*a24)
  if ((r[1]==2) & (r[2]==4) & (r[3]==2) & (r[4]==0))
    psi <-
      3*(2*a13^2*a22^2 + 16*a12*a13*a22*a23 + 4*a12^2*(2*a23^2 + a22*a33) + 
         a11*a22*(4*a23^2 + a22*a33))
  if ((r[1]==2) & (r[2]==4) & (r[3]==1) & (r[4]==1))
    psi <-
      3*(8*a12*a14*a22*a23 + 2*a13*a22*(a14*a22 + 4*a12*a24) + 
         4*a12^2*(2*a23*a24 + a22*a34) + a11*a22*(4*a23*a24 + a22*a34))
  if ((r[1]==2) & (r[2]==4) & (r[3]==0) & (r[4]==2))
    psi <-
      3*(2*a14^2*a22^2 + 16*a12*a14*a22*a24 + 4*a12^2*(2*a24^2 + a22*a44) + 
         a11*a22*(4*a24^2 + a22*a44))
  if ((r[1]==2) & (r[2]==3) & (r[3]==3) & (r[4]==0))
    psi <-
      3*(6*a13^2*a22*a23 + 6*a12*a13*(2*a23^2 + a22*a33) + 
         a23*(2*a11*a23^2 + 6*a12^2*a33 + 3*a11*a22*a33))
  if ((r[1]==2) & (r[2]==3) & (r[3]==2) & (r[4]==1))
    psi <-
      3*(2*a13^2*a22*a24 + 2*a12*a14*(2*a23^2 + a22*a33) + 
         4*a13*(a14*a22*a23 + 2*a12*a23*a24 + a12*a22*a34) + 
         2*a12^2*(a24*a33 + 2*a23*a34) + a11*(2*a23^2*a24 + a22*a24*a33 + 
                                              2*a22*a23*a34))
if ((r[1]==2) & (r[2]==3) & (r[3]==1) & (r[4]==2))
psi <-
3*(2*a14^2*a22*a23 + 4*a14*(a13*a22*a24 + 2*a12*a23*a24 + a12*a22*a34) + 
  2*a12*a13*(2*a24^2 + a22*a44) + 2*a12^2*(2*a24*a34 + a23*a44) + 
  a11*(2*a23*a24^2 + 2*a22*a24*a34 + a22*a23*a44))
if ((r[1]==2) & (r[2]==3) & (r[3]==0) & (r[4]==3))
psi <-
3*(6*a14^2*a22*a24 + 6*a12*a14*(2*a24^2 + a22*a44) + 
  a24*(2*a11*a24^2 + 6*a12^2*a44 + 3*a11*a22*a44))
if ((r[1]==2) & (r[2]==2) & (r[3]==4) & (r[4]==0))
psi <-
3*(16*a12*a13*a23*a33 + 4*a13^2*(2*a23^2 + a22*a33) + 
  a33*(4*a11*a23^2 + 2*a12^2*a33 + a11*a22*a33))
if ((r[1]==2) & (r[2]==2) & (r[3]==3) & (r[4]==1))
psi <-
3*(4*a12*a14*a23*a33 + 2*a12^2*a33*a34 + 2*a13^2*(2*a23*a24 + a22*a34) + 
  2*a13*(2*a14*a23^2 + a14*a22*a33 + 2*a12*a24*a33 + 4*a12*a23*a34) + 
  a11*(2*a23*a24*a33 + 2*a23^2*a34 + a22*a33*a34))
if ((r[1]==2) & (r[2]==2) & (r[3]==2) & (r[4]==2))
psi <-
2*a11*a24^2*a33 + 2*a14^2*(2*a23^2 + a22*a33) + 8*a11*a23*a24*a34 + 
 4*a12^2*a34^2 + 2*a11*a22*a34^2 + 
 8*a14*(2*a13*a23*a24 + a12*a24*a33 + a13*a22*a34 + 2*a12*a23*a34) + 
 2*a11*a23^2*a44 + 2*a12^2*a33*a44 + a11*a22*a33*a44 + 
 2*a13^2*(2*a24^2 + a22*a44) + 8*a12*a13*(2*a24*a34 + a23*a44)
if ((r[1]==2) & (r[2]==2) & (r[3]==1) & (r[4]==3))
psi <-
3*(2*a14^2*(2*a23*a24 + a22*a34) + 2*a12*(2*a13*a24 + a12*a34)*a44 + 
  2*a14*(2*a13*a24^2 + 4*a12*a24*a34 + a13*a22*a44 + 2*a12*a23*a44) + 
  a11*(2*a24^2*a34 + 2*a23*a24*a44 + a22*a34*a44))
  if ((r[1]==2) & (r[2]==2) & (r[3]==0) & (r[4]==4))
psi <-
3*(16*a12*a14*a24*a44 + 4*a14^2*(2*a24^2 + a22*a44) + 
  a44*(4*a11*a24^2 + 2*a12^2*a44 + a11*a22*a44))
if ((r[1]==2) & (r[2]==1) & (r[3]==5) & (r[4]==0))
psi <-
15*a33*(4*a13^2*a23 + 2*a12*a13*a33 + a11*a23*a33)
if ((r[1]==2) & (r[2]==1) & (r[3]==4) & (r[4]==1))
psi <-
3*(8*a13*a33*(a14*a23 + a12*a34) + 4*a13^2*(a24*a33 + 2*a23*a34) + 
  a33*(2*a12*a14*a33 + a11*a24*a33 + 4*a11*a23*a34))
if ((r[1]==2) & (r[2]==1) & (r[3]==3) & (r[4]==2))
psi <-
3*(2*a14^2*a23*a33 + 4*a14*(a13*a24*a33 + 2*a13*a23*a34 + a12*a33*a34) + 
  2*a13^2*(2*a24*a34 + a23*a44) + 2*a12*a13*(2*a34^2 + a33*a44) + 
  a11*(2*a24*a33*a34 + 2*a23*a34^2 + a23*a33*a44))
if ((r[1]==2) & (r[2]==1) & (r[3]==2) & (r[4]==3))
psi <-
3*(2*a14^2*(a24*a33 + 2*a23*a34) + 2*a13*(a13*a24 + 2*a12*a34)*a44 + 
  2*a14*(4*a13*a24*a34 + 2*a12*a34^2 + 2*a13*a23*a44 + a12*a33*a44) + 
  a11*(2*a24*a34^2 + a24*a33*a44 + 2*a23*a34*a44))