📄 bivar.f90
字号:
p21 = 0.0E+00
p23 = -zuu(2)+zuu(1)
p22 = -1.5E+00*p23
itpv = it0
!
! Convert XII and YII to UV system.
!
50 continue
dx = xii-x0
dy = yii-y0
u = ap*dx+bp*dy
v = cp*dx+dp*dy
!
! Evaluate the polynomial.
!
p0 = p00+v*(p01+v*(p02+v*(p03+v*(p04+v*p05))))
p1 = p10+v*(p11+v*(p12+v*p13))
p2 = p20+v*(p21+v*(p22+v*p23))
zii = p0+u*(p1+u*p2)
return
!
! Calculation of ZII by extrapolation in the triangle.
! Check if the necessary coefficients have been calculated.
!
60 continue
if ( it0 /= itpv ) then
!
! Load coordinate and partial derivative values at the vertex of the triangle.
!
jipl = 3*il2-2
idp = ipl(jipl)
x0 = xd(idp)
y0 = yd(idp)
z0 = zd(idp)
jpdd = 5*(idp-1)
do kpd = 1, 5
jpdd = jpdd+1
pd(kpd) = pdd(jpdd)
end do
!
! Calculate the coefficients of the polynomial.
!
p00 = z0
p10 = pd(1)
p01 = pd(2)
p20 = 0.5E+00*pd(3)
p11 = pd(4)
p02 = 0.5E+00*pd(5)
itpv = it0
end if
!
! Convert XII and YII to UV system.
!
u = xii-x0
v = yii-y0
!
! Evaluate the polynomial.
!
p0 = p00+v*(p01+v*p02)
p1 = p10+v*p11
zii = p0+u*(p1+u*p20)
return
end
subroutine idsfft ( md, ndp, xd, yd, zd, nxi, nyi, nzi, xi, yi, zi, iwk, wk )
!
!*******************************************************************************
!
!! IDSFFT fits a smooth surface Z(X,Y) given irregular (X,Y,Z) data.
!
!
! Discussion:
!
! IDSFFT performs smooth surface fitting when the projections of the
! data points in the (X,Y) plane are irregularly distributed.
!
! Special conditions:
!
! Inadequate work space IWK and WK may may cause incorrect results.
!
! The data points must be distinct and their projections in the XY
! plane must not be collinear, otherwise an error return occurs.
!
! Parameters:
!
! Input, integer MD, mode of computation (must be 1, 2, or 3,
! else an error return will occur).
!
! 1, if this is the first call to this routine, or if the value of
! NDP has been changed from the previous call, or if the contents of
! the XD or YD arrays have been changed from the previous call.
!
! 2, if the values of NDP and the XD, YD arrays are unchanged from
! the previous call, but new values for XI, YI are being used. If
! MD = 2 and NDP has been changed since the previous call to IDSFFT,
! an error return occurs.
!
! 3, if the values of NDP, NXI, NYI, XD, YD, XI, YI are unchanged
! from the previous call, i.e. if the only change on input to idsfft
! is in the ZD array. If MD = 3 and NDP, nxi or nyi has been changed
! since the previous call to idsfft, an error return occurs.
!
! Between the call with MD = 2 or MD = 3 and the preceding call, the
! iwk and wk work arrays should not be disturbed.
!
! Input, integer NDP, the number of data points. NDP must be at least 4.
!
! Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data.
!
! Input, real ZD(NDP), the data values.
!
! Input, integer NXI, NYI, the number of output grid points in the
! X and Y directions. NXI and NYI must each be at least 1.
!
! Input, integer NZI, the first dimension of ZI. NZI must be at
! least NXI.
!
! Input, real XI(NXI), YI(NYI), the X and Y coordinates of the grid
! points.
!
! Workspace, integer IWK(31*NDP+NXI*NYI).
!
! Workspace, real WK(6*NDP).
!
! Output, real ZI(NZI,NYI), contains the interpolated Z values at the
! grid points.
!
implicit none
!
integer ndp
integer nxi
integer nyi
integer nzi
!
real ap
real bp
real cp
real dp
integer il1
integer il2
integer iti
integer itpv
integer iwk(31*ndp + nxi*nyi)
integer ixi
integer iyi
integer izi
integer jig0mn
integer jig0mx
integer jig1mn
integer jig1mx
integer jigp
integer jngp
integer jwigp
integer jwigp0
integer jwipl
integer jwipt
integer jwiwl
integer jwiwp
integer jwngp
integer jwngp0
integer jwwpd
integer md
integer ngp0
integer ngp1
integer nl
integer nngp
integer nt
real p00
real p01
real p02
real p03
real p04
real p05
real p10
real p11
real p12
real p13
real p14
real p20
real p21
real p22
real p23
real p30
real p31
real p32
real p40
real p41
real p50
real wk(6*ndp)
real x0
real xd(ndp)
real xi(nxi)
real y0
real yd(ndp)
real yi(nyi)
real zd(ndp)
real zi(nzi,nyi)
!
save /idpt/
!
common /idpt/ itpv,x0,y0,ap,bp,cp,dp, &
p00,p10,p20,p30,p40,p50,p01,p11,p21,p31,p41, &
p02,p12,p22,p32,p03,p13,p23,p04,p14,p05
!
! Error check.
!
if ( md < 1 .or. md > 3 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write(*,*)' Input parameter MD out of range.'
stop
end if
if ( ndp < 4 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write ( *, '(a)' ) ' Input parameter NDP out of range.'
stop
end if
if ( nxi < 1 .or. nyi < 1 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write ( *, '(a)' ) ' Input parameter NXI or NYI out of range.'
stop
end if
if ( nxi > nzi ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write ( *, '(a)' ) ' Input parameter NZI is less than NXI.'
stop
end if
if ( md <= 1 ) then
iwk(1) = ndp
else
if ( ndp /= iwk(1) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write ( *, '(a)' ) ' MD = 2 or 3 but ndp was changed since last call.'
stop
end if
end if
if ( md <= 2 ) then
iwk(3) = nxi
iwk(4) = nyi
else
if ( nxi /= iwk(3) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write ( *, '(a)' ) 'MD = 3 but nxi was changed since last call.'
stop
end if
if ( nyi /= iwk(4) ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDSFFT - Fatal error!'
write ( *, '(a)' ) ' MD = 3 but nyi was changed since last call.'
stop
end if
end if
!
! Allocation of storage areas in the IWK array.
!
jwipt = 16
jwiwl = 6*ndp+1
jwngp0 = jwiwl-1
jwipl = 24*ndp+1
jwiwp = 30*ndp+1
jwigp0 = 31*ndp
jwwpd = 5*ndp+1
!
! Triangulate the XY plane.
!
if ( md == 1 ) then
call idtang ( ndp, xd, yd, nt, iwk(jwipt), nl, iwk(jwipl), &
iwk(jwiwl), iwk(jwiwp), wk )
iwk(5) = nt
iwk(6) = nl
if ( nt == 0 ) then
return
end if
else
nt = iwk(5)
nl = iwk(6)
end if
!
! Sort output grid points in ascending order of the triangle
! number and the border line segment number.
!
if ( md <= 2 ) then
call idgrid ( xd, yd, nt, iwk(jwipt), nl, iwk(jwipl), nxi, &
nyi, xi, yi, iwk(jwngp0+1), iwk(jwigp0+1) )
end if
!
! Estimate partial derivatives at all data points.
!
call idpdrv ( ndp, xd, yd, zd, nt, iwk(jwipt), wk, wk(jwwpd) )
!
! Interpolate the ZI values.
!
itpv = 0
jig0mx = 0
jig1mn = nxi*nyi+1
nngp = nt+2*nl
do jngp = 1, nngp
iti = jngp
if ( jngp > nt ) then
il1 = (jngp-nt+1)/2
il2 = (jngp-nt+2)/2
if(il2>nl) then
il2 = 1
end if
iti = il1*(nt+nl)+il2
end if
jwngp = jwngp0+jngp
ngp0 = iwk(jwngp)
if ( ngp0 /= 0 ) then
jig0mn = jig0mx+1
jig0mx = jig0mx+ngp0
do jigp = jig0mn, jig0mx
jwigp = jwigp0+jigp
izi = iwk(jwigp)
iyi = (izi-1)/nxi+1
ixi = izi-nxi*(iyi-1)
call idptip ( ndp, xd, yd, zd, nt, iwk(jwipt), nl, iwk(jwipl), &
wk, iti, xi(ixi), yi(iyi), zi(ixi,iyi) )
end do
end if
jwngp = jwngp0+2*nngp+1-jngp
ngp1 = iwk(jwngp)
if ( ngp1 /= 0 ) then
jig1mx = jig1mn-1
jig1mn = jig1mn-ngp1
do jigp = jig1mn, jig1mx
jwigp = jwigp0+jigp
izi = iwk(jwigp)
iyi = (izi-1)/nxi+1
ixi = izi-nxi*(iyi-1)
call idptip ( ndp, xd, yd, zd, nt, iwk(jwipt), nl, iwk(jwipl), &
wk, iti, xi(ixi), yi(iyi), zi(ixi,iyi) )
end do
end if
end do
return
end
subroutine idtang ( ndp, xd, yd, nt, ipt, nl, ipl, iwl, iwp, wk )
!
!*******************************************************************************
!
!! IDTANG performs triangulation.
!
!
! Discussion:
!
! The routine divides the XY plane into a number of triangles according to
! given data points in the plane, determines line segments that form
! the border of data area, and determines the triangle numbers
! corresponding to the border line segments.
!
! At completion, point numbers of the vertexes of each triangle
! are listed counter-clockwise. Point numbers of the end points
! of each border line segment are listed counter-clockwise,
! listing order of the line segments being counter-clockwise.
!
! Parameters:
!
! Input, integer NDP, the number of data points.
!
! Input, real XD(NDP), YD(NDP), the X and Y coordinates of the data.
!
! Output, integer NT, the number of triangles,
!
! Output, integer IPT(6*NDP-15), where the point numbers of the
! vertexes of the IT-th triangle are to be stored as entries
! 3*IT-2, 3*IT-1, and 3*IT, for IT = 1 to NT.
!
! Output, integer NL, the number of border line segments.
!
! Output, integer IPL(6*NDP), where the point numbers of the end
! points of the (il)th border line segment and its respective triangle
! number are to be stored as the (3*il-2)nd, (3*il-1)st, and (3*il)th
! elements, il = 1,2,..., nl.
!
! Workspace, integer IWL(18*NDP),
!
! Workspace, integer IWP(NDP),
!
! Workspace, real WK(NDP).
!
implicit none
!
integer ndp
!
real dsqf
real dsqi
real dsqmn
real, parameter :: epsln = 1.0E-06
integer idxchg
integer il
integer ilf
integer iliv
integer ilt3
integer ilvs
integer ip
integer ip1
integer ip1p1
integer ip2
integer ip3
integer ipl(6*ndp)
integer ipl1
integer ipl2
integer iplj1
integer iplj2
integer ipmn1
integer ipmn2
integer ipt(6*ndp-15)
integer ipt1
integer ipt2
integer ipt3
integer ipti
integer ipti1
integer ipti2
integer irep
integer it
integer it1t3
integer it2t3
integer itf(2)
integer its
integer itt3
integer itt3r
integer iwl(18*ndp)
integer iwp(ndp)
integer ixvs
integer ixvspv
integer jl1
integer jl2
integer jlt3
integer jp
integer jp1
integer jp2
integer jpc
integer jpmn
integer jpmx
integer jwl
integer jwl1
integer jwl1mn
integer nl
integer nl0
integer nlf
integer nlfc
integer nlft2
integer nln
integer nlnt3
integer nlsh
integer nlsht3
integer nlt3
integer, parameter :: nrep = 100
integer nt
integer nt0
integer ntf
integer ntt3
integer ntt3p3
real sp
real spdt
real u1
real u2
real u3
real v1
real v2
real v3
real vp
real vpdt
real wk(ndp)
real x1
real x2
real x3
real xd(ndp)
real xdmp
real y1
real y2
real y3
real yd(ndp)
real ydmp
!
! Statement functions
!
dsqf(u1,v1,u2,v2) = (u2-u1)**2+(v2-v1)**2
spdt(u1,v1,u2,v2,u3,v3) = (u2-u1)*(u3-u1)+(v2-v1)*(v3-v1)
vpdt(u1,v1,u2,v2,u3,v3) = (v3-v1)*(u2-u1)-(u3-u1)*(v2-v1)
!
! Preliminary processing
!
if ( ndp < 4 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDTANG - Fatal error!'
write ( *, '(a)' ) ' Input parameter NDP out of range.'
stop
end if
!
! Determine IPMN1 and IPMN2, the closest pair of data points.
!
dsqmn = dsqf(xd(1),yd(1),xd(2),yd(2))
ipmn1 = 1
ipmn2 = 2
do ip1 = 1, ndp-1
x1 = xd(ip1)
y1 = yd(ip1)
ip1p1 = ip1+1
do ip2 = ip1p1, ndp
dsqi = dsqf(x1,y1,xd(ip2),yd(ip2))
if ( dsqi == 0.0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'IDTANG - Fatal error!'
write ( *, '(a)' ) ' Two of the input data points are identical.'
stop
end if
if(dsqi<dsqmn) then
dsqmn = dsqi
ipmn1 = ip1
ipmn2 = ip2
end if
end do
end do
!
! Compute the midpoint of the closest two data points.
!
xdmp = (xd(ipmn1)+xd(ipmn2)) / 2.0E+00
ydmp = (yd(ipmn1)+yd(ipmn2)) / 2.0E+00
!
! Sort the other (NDP-2) data points in ascending order of
! distance from the midpoint and store the sorted data point
! numbers in the IWP array.
!
jp1 = 2
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -