V3FIT
dnherm3.f
1  subroutine dnherm3(x,nx,y,ny,z,nz,fherm,nf2,nf3,
2  > ilinx,iliny,ilinz,ier)
3 C
4 C create a data set for Hermite interpolation, based on simple
5 C numerical differentiation using the given grid. The grid should
6 C be "approximately" evenly spaced.
7 C
8 C 3d routine
9 C
10 C input:
11 C
12 C nf2.get.nx and nf3.ge.ny expected!
13 C
14  integer nx,ny,nz,nf2,nf3 ! array dimensions
15  real x(nx) ! x coordinate array
16  real y(ny) ! y coordinate array
17  real z(nz) ! z coordinate array
18  real fherm(0:7,nf2,nf3,nz) ! data/Hermite array
19 C
20 C fherm(0,i,j,k) = function value f at x(i),y(j),z(k) **on input**
21 C
22 C fherm(1,i,j,k) = derivative df/dx at x(i),y(j),z(k) **on output**
23 C fherm(2,i,j,k) = derivative df/dy at x(i),y(j),z(k) **on output**
24 C fherm(3,i,j,k) = derivative df/dz at x(i),y(j),z(k) **on output**
25 C fherm(4,i,j,k) = derivative d2f/dxdy at x(i),y(j),z(k) **on output**
26 C fherm(5,i,j,k) = derivative d2f/dxdz at x(i),y(j),z(k) **on output**
27 C fherm(6,i,j,k) = derivative d2f/dydz at x(i),y(j),z(k) **on output**
28 C fherm(7,i,j,k) = derivative d3f/dxdydz at x(i),y(j),z(k) **on output**
29 C
30 C additional outputs:
31 C
32 C ilinx = 1 if x is evenly spaced (approximately)
33 C iliny = 1 if y is evenly spaced (approximately)
34 C ilinz = 1 if z is evenly spaced (approximately)
35 C
36 C ier = 0 if there are no errors
37 C
38 C note possible errors: x y or z NOT strict ascending
39 C nf2.lt.nx .or. nf3.lt.ny
40 C
41 C----------------------------
42 C
43 C
44 C error check
45 C
46  ier=0
47 c
48  call splinck(x,nx,ilinx,1.0e-3,ierx)
49  if(ierx.ne.0) ier=2
50 c
51  if(ier.eq.2) then
52  write(6,'('' ?dnherm3: x axis not strict ascending'')')
53  endif
54 c
55  call splinck(y,ny,iliny,1.0e-3,iery)
56  if(iery.ne.0) ier=3
57 c
58  if(ier.eq.3) then
59  write(6,'('' ?dnherm3: y axis not strict ascending'')')
60  endif
61 c
62  call splinck(z,nz,ilinz,1.0e-3,ierz)
63  if(ierz.ne.0) ier=4
64 c
65  if(ier.eq.4) then
66  write(6,'('' ?dnherm3: z axis not strict ascending'')')
67  endif
68 c
69  if(nf2.lt.nx) then
70  ier=5
71  write(6,*) '?dnherm3: fherm (x) array dimension too small.'
72  endif
73 c
74  if(nf3.lt.ny) then
75  ier=6
76  write(6,*) '?dnherm3: fherm (y) array dimension too small.'
77  endif
78 C
79  if(ier.ne.0) return
80 C
81 C error check OK
82 C
83  do iz=1,nz
84 c
85  izp=min(nz,iz+1)
86  izm=max(1,iz-1)
87 c
88  do iy=1,ny
89 c
90  iyp=min(ny,iy+1)
91  iym=max(1,iy-1)
92 c
93  do ix=1,nx
94 c
95  ixp=min(nx,ix+1)
96  ixm=max(1,ix-1)
97 c
98 c x diffs in vicinity
99 c
100  zd=(fherm(0,ixp,iy,iz)-fherm(0,ixm,iy,iz))/
101  > (x(ixp)-x(ixm))
102 c
103  fherm(1,ix,iy,iz)=zd
104 c
105 c y diffs in vicinity
106 c
107  zd=(fherm(0,ix,iyp,iz)-fherm(0,ix,iym,iz))/
108  > (y(iyp)-y(iym))
109 c
110  fherm(2,ix,iy,iz)=zd
111 c
112 c z diffs in vicinity
113 c
114  zd=(fherm(0,ix,iy,izp)-fherm(0,ix,iy,izm))/
115  > (z(izp)-z(izm))
116 c
117  fherm(3,ix,iy,iz)=zd
118 c
119 c xy cross derivs
120 c
121  fherm(4,ix,iy,iz)=
122  > (fherm(0,ixp,iyp,iz)-fherm(0,ixm,iyp,iz)
123  > -fherm(0,ixp,iym,iz)+fherm(0,ixm,iym,iz))/
124  > ((x(ixp)-x(ixm))*(y(iyp)-y(iym)))
125 c
126 c xz cross derivs
127 c
128  fherm(5,ix,iy,iz)=
129  > (fherm(0,ixp,iy,izp)-fherm(0,ixm,iy,izp)
130  > -fherm(0,ixp,iy,izm)+fherm(0,ixm,iy,izm))/
131  > ((x(ixp)-x(ixm))*(z(izp)-z(izm)))
132 c
133 c yz cross derivs
134 c
135  fherm(6,ix,iy,iz)=
136  > (fherm(0,ix,iyp,izp)-fherm(0,ix,iym,izp)
137  > -fherm(0,ix,iyp,izm)+fherm(0,ix,iym,izm))/
138  > ((y(iyp)-y(iym))*(z(izp)-z(izm)))
139 c
140 c xyz cross deriv
141 c
142  fherm(7,ix,iy,iz)=
143  > ((fherm(0,ixp,iyp,izp)-fherm(0,ixp,iym,izp)
144  > -fherm(0,ixp,iyp,izm)+fherm(0,ixp,iym,izm))-
145  > (fherm(0,ixm,iyp,izp)-fherm(0,ixm,iym,izp)
146  > -fherm(0,ixm,iyp,izm)+fherm(0,ixm,iym,izm)))/
147  > ((x(ixp)-x(ixm))*(y(iyp)-y(iym))*(z(izp)-z(izm)))
148 c
149  enddo
150  enddo
151  enddo
152 C
153  return
154  end