! test_4d.f90  just the test function and derivatives
!              fourth order and through fourth derivative
! u(x,y,z,t) and ux(x,y,z,t) etc.  many of them

! every possible term with a unique constant coefficient

module test_4d
contains

function u(x, y, z, t) result(val)
  implicit none
  double precision :: x, y, z, t, val
  val = x*x*x*x +  2.0*x*x*x*y +  3.0*x*x*x*z +  4.0*x*x*x*t +  &
     5.0*x*x*y*y +  6.0*x*x*y*z +  7.0*x*x*y*t +  8.0*x*x*z*z + &
     9.0*x*x*z*t + 10.0*x*x*t*t + 11.0*x*y*y*y + 12.0*x*y*y*z + &
    13.0*x*y*y*t + 14.0*x*y*z*z + 15.0*x*y*z*t + 16.0*x*y*t*t + &
    17.0*x*z*z*z + 18.0*x*z*z*t + 19.0*x*z*t*t + 20.0*x*t*t*t + &
    21.0*y*y*y*y + 22.0*y*y*y*z + 23.0*y*y*y*t + 24.0*y*y*z*z + &
    25.0*y*y*z*t + 26.0*y*y*t*t + 27.0*y*z*z*z + 28.0*y*z*z*t + &
    29.0*y*z*t*t + 30.0*y*t*t*t + 31.0*z*z*z*z + 32.0*z*z*z*t + &
    33.0*z*z*t*t + 34.0*z*t*t*t + 35.0*t*t*t*t +                &
    36.0*x*x*x + 37.0*x*x*y + 38.0*x*x*z + 39.0*x*x*t + 40.0*x*y*y + &
    41.0*x*y*z + 42.0*x*y*t + 43.0*x*z*z + 44.0*x*z*t + 45.0*x*t*t + &
    46.0*y*y*y + 47.0*y*y*z + 48.0*y*y*t + 49.0*y*z*z + 50.0*y*z*t + &
    51.0*y*t*t + 52.0*z*z*z + 53.0*z*z*t + 54.0*z*t*t + 55.0*t*t*t + &
    56.0*x*x + 57.0*x*y + 58.0*x*z + 59.0*x*t + 60.0*y*y +      &
    61.0*y*z + 62.0*y*t + 63.0*z*z + 64.0*z*t + 65.0*t*t +      &
    66.0*x + 67.0*y + 68.0*z + 69.0*t + 70.0
end function u 

! now, for testing, the 69 possible non zero derivatives

function Uxxxx(x, y, z, t) result(val)
  implicit none
  double precision :: x, y, z, t, val
  val = 24.0;
end function Uxxxx

function Uxxxy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 12.0;
end function Uxxxy;

function Uxxxz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 18.0;
end function Uxxxz;

function Uxxxt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 24.0;
end function Uxxxt;

function Uxxyy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 20.0;
end function Uxxyy;

function Uxxyz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 12.0;
end function Uxxyz;

function Uxxyt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 14.0;
end function Uxxyt;

function Uxxzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 32.0;
end function Uxxzz;

function Uxxzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 18.0;
end function Uxxzt;

function Uxxtt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 40.0;
end function Uxxtt;

function Uxyyy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 66.0;
end function Uxyyy;

function Uxyyz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 24.0;
end function Uxyyz;

function Uxyyt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 26.0;
end function Uxyyt;

function Uxyzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 28.0;
end function Uxyzz;

function Uxyzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 15.0;
end function Uxyzt;

function Uxytt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 32.0;
end function Uxytt;

function Uxzzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 102.0;
end function Uxzzz;

function Uxzzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 36.0;
end function Uxzzt;

function Uxztt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 38.0;
end function Uxztt;

function Uxttt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 120.0;
end function Uxttt;

function Uyyyy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 504.0;
end function Uyyyy;

function Uyyyz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 132.0;
end function Uyyyz;

function Uyyyt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 138.0;
end function Uyyyt;

function Uyyzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 96.0;
end function Uyyzz;

function Uyyzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 50.0;
end function Uyyzt;

function Uyytt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 104.0;
end function Uyytt;

function Uyzzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 162.0;
end function Uyzzz;

function Uyzzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 56.0;
end function Uyzzt;

function Uyztt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 58.0;
end function Uyztt;

function Uyttt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 180.0;
end function Uyttt;

function Uzzzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 744.0;
end function Uzzzz;

function Uzzzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 192.0;
end function Uzzzt;

function Uzztt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 132.0;
end function Uzztt;

function Uzttt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 204.0;
end function Uzttt;

function Utttt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 840.0;
end function Utttt;

function Uxxx(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 12.0*y+18.0*z+216.0+24.0*x+24.0*t;
end function Uxxx;

function Uxxy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 12.0*x+14.0*t+74.0+20.0*y+12.0*z;
end function Uxxy;

function Uxxz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 18.0*t+18.0*x+76.0+32.0*z+12.0*y;
end function Uxxz;

function Uxxt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 18.0*z+78.0+14.0*y+24.0*x+40.0*t;
end function Uxxt;

function Uxyy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 26.0*t+80.0+20.0*x+24.0*z+66.0*y;
end function Uxyy;

function Uxyz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 41.0+28.0*z+15.0*t+24.0*y+12.0*x;
end function Uxyz;

function Uxyt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 26.0*y+42.0+14.0*x+32.0*t+15.0*z;
end function Uxyt;

function Uxzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 86.0+36.0*t+28.0*y+102.0*z+32.0*x;
end function Uxzz;

function Uxzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 38.0*t+18.0*x+36.0*z+44.0+15.0*y;
end function Uxzt;

function Uxtt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 38.0*z+120.0*t+90.0+32.0*y+40.0*x;
end function Uxtt;

function Uyyy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 138.0*t + 276.0 + 66.0*x + 132.0*z + 504.0*y;
end function Uyyy;

function Uyyz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 94.0+96.0*z+50.0*t+132.0*y+24.0*x;
end function Uyyz;

function Uyyt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 138.0*y+96.0+26.0*x+104.0*t+50.0*z;
end function Uyyt;

function Uyzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 98.0+56.0*t+96.0*y+162.0*z+28.0*x;
end function Uyzz;

function Uyzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 58.0*t+15.0*x+56.0*z+50.0+50.0*y;
end function Uyzt;

function Uytt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 58.0*z+180.0*t+102.0+104.0*y+32.0*x;
end function Uytt;

function Uzzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 312.0+192.0*t+162.0*y+744.0*z+102.0*x;
end function Uzzz;

function Uzzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 132.0*t+36.0*x+192.0*z+106.0+56.0*y;
end function Uzzt;

function Uztt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 132.0*z+204.0*t+108.0+58.0*y+38.0*x;
end function Uztt;

function Uttt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 204.0*z+840.0*t+330.0+180.0*y+120.0*x;
end function Uttt;

function Uxx(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 24.0*x*t+14.0*y*t+12.0*y*z+20.0*t*t+78.0*t+12.0*x*y+18.0*x*z+76.0*z+&
        74.0*y+216.0*x+10.0*y*y+16.0*z*z+112.0+18.0*z*t+12.0*x*x;
end function Uxx;

function Uxy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 80.0*y+12.0*x*z+33.0*y*y+6.0*x*x+74.0*x+15.0*z*t+16.0*t*t+24.0*y*z+&
        26.0*y*t+41.0*z+57.0+42.0*t+20.0*x*y+14.0*z*z+14.0*x*t;
end function Uxy;

function Uxz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 51.0*z*z+12.0*y*y+12.0*x*y+58.0+15.0*y*t+18.0*x*t+86.0*z+44.0*t+&
        76.0*x+9.0*x*x+19.0*t*t+41.0*y+28.0*y*z+36.0*z*t+32.0*x*z;
end function Uxz;

function Uxt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 14.0*x*y+18.0*x*z+40.0*x*t+78.0*x+44.0*z+12.0*x*x+90.0*t+59.0+&
        38.0*z*t+15.0*y*z+42.0*y+18.0*z*z+60.0*t*t+13.0*y*y+32.0*y*t;
end function Uxt;

function Uyy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 26.0*x*t+276.0*y+52.0*t*t+50.0*z*t+24.0*x*z+138.0*y*t+10.0*x*x+&
        80.0*x+48.0*z*z+94.0*z+120.0+96.0*t+252.0*y*y+66.0*x*y+132.0*y*z;
end function Uyy;

function Uyz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 56.0*z*t+24.0*x*y+15.0*x*t+61.0+98.0*z+50.0*t+41.0*x+6.0*x*x+&
        81.0*z*z+66.0*y*y+29.0*t*t+28.0*x*z+94.0*y+50.0*y*t+96.0*y*z;
end function Uyz;

function Uyt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 26.0*x*y+32.0*x*t+42.0*x+15.0*x*z+50.0*z+102.0*t+62.0+58.0*z*t+&
        28.0*z*z+104.0*y*t+7.0*x*x+96.0*y+50.0*y*z+69.0*y*y+90.0*t*t;
end function Uyt;

function Uzz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 126.0+28.0*x*y+102.0*x*z+312.0*z+106.0*t+36.0*x*t+66.0*t*t+86.0*x+&
        192.0*z*t+16.0*x*x+56.0*y*t+98.0*y+162.0*y*z+372.0*z*z+48.0*y*y;
end function Uzz;

function Uzt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 15.0*x*y+44.0*x+58.0*y*t+106.0*z+9.0*x*x+36.0*x*z+132.0*z*t+108.0*t+&
        64.0+102.0*t*t+96.0*z*z+38.0*x*t+56.0*y*z+50.0*y+25.0*y*y;
end function Uzt;

function Utt(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val = 32.0*x*y+90.0*x+180.0*y*t+108.0*z+20.0*x*x+38.0*x*z+204.0*z*t+&
        330.0*t+130.0+420.0*t*t+66.0*z*z+120.0*x*t+58.0*y*z+102.0*y+52.0*y*y;
end function Utt;

function Ux(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val =  40.0*y*y+14.0*x*y*t+58.0*z+12.0*x*y*z+10.0*x*y*y+16.0*x*z*z+       &
         78.0*x*t+43.0*z*z+16.0*y*t*t+44.0*z*t+12.0*x*x*t+18.0*x*z*t+       &
         19.0*z*t*t+4.0*x*x*x+76.0*x*z+74.0*x*y+45.0*t*t+59.0*t+6.0*x*x*y+  &
         11.0*y*y*y+20.0*t*t*t+18.0*z*z*t+9.0*x*x*z+20.0*x*t*t+15.0*y*z*t+  &
         108.0*x*x+112.0*x+41.0*y*z+14.0*y*z*z+57.0*y+17.0*z*z*z+12.0*y*y*z+&
         42.0*y*t+13.0*y*y*t+66.0;
end function Ux;

function Uy(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val =  138.0*y*y+26.0*x*y*t+61.0*z+24.0*x*y*z+33.0*x*y*y+14.0*x*z*z+     &
         42.0*x*t+49.0*z*z+52.0*y*t*t+50.0*z*t+7.0*x*x*t+15.0*x*z*t+       &
         29.0*z*t*t+2.0*x*x*x+41.0*x*z+80.0*x*y+51.0*t*t+62.0*t+10.0*x*x*y+&
         84.0*y*y*y+30.0*t*t*t+28.0*z*z*t+6.0*x*x*z+16.0*x*t*t+50.0*y*z*t+ &
         37.0*x*x+57.0*x+94.0*y*z+48.0*y*z*z+120.0*y+27.0*z*z*z+66.0*y*y*z+&
         96.0*y*t+67.0+69.0*y*y*t;
end function Uy;

function Uz(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val =  47.0*y*y+15.0*x*y*t+126.0*z+28.0*x*y*z+12.0*x*y*y+51.0*x*z*z+     &
         44.0*x*t+156.0*z*z+29.0*y*t*t+106.0*z*t+9.0*x*x*t+36.0*x*z*t+     &
         66.0*z*t*t+3.0*x*x*x+86.0*x*z+41.0*x*y+54.0*t*t+64.0*t+6.0*x*x*y+ &
         22.0*y*y*y+34.0*t*t*t+96.0*z*z*t+16.0*x*x*z+19.0*x*t*t+56.0*y*z*t+&
         38.0*x*x+58.0*x+98.0*y*z+81.0*y*z*z+61.0*y+124.0*z*z*z+48.0*y*y*z+&
         50.0*y*t+68.0+25.0*y*y*t;
end function Uz;

function Ut(x, y, z, t) result(val)
  double precision :: x, y, z, t, val
  val =  48.0*y*y+32.0*x*y*t+64.0*z+15.0*x*y*z+13.0*x*y*y+18.0*x*z*z+     &
         90.0*x*t+53.0*z*z+90.0*y*t*t+108.0*z*t+20.0*x*x*t+38.0*x*z*t+    &
         102.0*z*t*t+4.0*x*x*x+44.0*x*z+42.0*x*y+165.0*t*t+130.0*t+       &
         7.0*x*x*y+23.0*y*y*y+140.0*t*t*t+66.0*z*z*t+9.0*x*x*z+60.0*x*t*t+&
         58.0*y*z*t+39.0*x*x+69.0+59.0*x+50.0*y*z+28.0*y*z*z+62.0*y+      &
         32.0*z*z*z+25.0*y*y*z+102.0*y*t+52.0*y*y*t;
end function Ut;

end module Test_4d;