[原创]抛砖引玉(关于二维插值和坐标转换)，请高手指点

bosshawk

问题：
对于一个网格，已知四顶点的坐标和四顶点处的流体密度值，求网格内某点(给定坐标)的密度值。
如图，四顶点标记为1，2，3，4，网格内点标记为P。注意1，2，3，4点可以组成任意凸四边形。
        (1)
        /\
       /  \
      /    \
     /  。  \
(2) \   (P)/ (3)
      \    /
       \  /
        \/
        (4)
设想一个逻辑坐标系，在该坐标系内四顶点的坐标分别为：1(1，1)，2(0，1)，3(1，0)，4(0，0)
如果能够根据实际坐标值求得P点在逻辑坐标系内的坐标P(a，b) (|a|<=1,|b|<=1)，那么可以采用双线性插值法[Rho(x) 表示 x 点的密度]：
Rho(P) = (1-a)(1-b) x Rho(4) + (1-a)b x Rho(2) + a(1-b) x Rho(3) + ab x Rho(1)
现在的问题是：已知实际坐标(x1,y1),(x2,y2),(x3,y3),(x4,y4),(xP,yP),怎样求a,b ？
下面是一段FORTRAN格式的程序，用于求a,b
      subroutine Calcab()
      d_max      =  1.d-8  !收敛条件
      maxiter    =  20     !最大循环次数
      a          =  0.5
      b          =  0.5
c     a,b可以给定任意0到1之间的初值

      d_x_xie    = 0.5*(x1-x2+x3-x4)
      d_y_xie    = 0.5*(y1-y2+y3-y4)
      d_x_eta    = 0.5*(x1-x3+x2-x4)
      d_y_eta    = 0.5*(y1-y3+y2-y4)

      deta  = d_x_xie*d_y_eta - d_y_xie*d_x_eta
      rden  = 1.0/deta
      test  = d_max*abs(deta)
      iter  = 0
   11 continue

      a1    = 1.0 - a
      b1    = 1.0 - b
      s1    = a*b
      s2    = a1*b
      s3    = a*b1
      s4    = a1*b1

      d_x   = x1*s1 + x2*s2 + x3*s3 + x4*s4 - xp
      d_y   = y1*s1 + y2*s2 + y3*s3 + y4*s4 - yp

      a = a + (d_y*d_x_eta - d_x*d_y_eta)*rden
      b = b + (d_x*d_y_xie - d_y*d_x_xie)*rden

      if (d_x**2 + d_y**2 .le. test) go to 10
      iter  = iter + 1
      if (iter.lt.maxiter) goto 11
   10 continue

      write(*,*) 'a:', a, '   b:',b

      if (iter.ge.maxiter) then
        write(*,*) '未收敛'
        stop
      endif
      end
请指点：这段程序是我从其他代码里抠出来的，只知道作用，不知道原理，请哪一位知道的指点一二，谢谢！

boling2000

三线性插值我也自己做了并且一直在用，好像没有什么问题呀？
——抱歉没有时间帮你查看程序，其实是一个非线性迭代，应该很快就收敛了。

bosshawk

对不起，没有说清楚，我关心的是a,b是怎么求出来的，也就是坐标怎么变换的。
那段程序没有问题，我想知道的是它依据的原理是什么，能否指引一下该看什么书或者资料？ 谢谢！！