Return to interp.f CVS log

Up to [Research Unix] / researchv10no / cmd / view2d / num

Annotation of researchv10no/cmd/view2d/num/interp.f, revision 1.1.1.1

1.1       root        1:       subroutine interp(u,d1,d2,d3)
                      2: c     interpolation in a triangle
                      3: c     input:  u = point to be interpolated
                      4: c             d1..3 = 3 points already known
                      5: c                 first subscript: xyz, second: point index
                      6: c     (output is in third component of u)
                      7:       real u(3)
                      8:       real d1(3),d2(3),d3(3)
                      9:       integer i,ipvt(3),info
                     10:       real a(3,3), r(3)
                     11:       do 100 i=1,3
                     12:         a(i,1)=1
                     13:  100   continue
                     14:       do 101 i=2,3
                     15:         a(1,i)=d1(i-1)
                     16:         a(2,i)=d2(i-1)
                     17:         a(3,i)=d3(i-1)
                     18:  101   continue
                     19:       r(1)=d1(3)
                     20:       r(2)=d2(3)
                     21:       r(3)=d3(3)
                     22:       call sgefa(a,3,3,ipvt,info)
                     23:       if(info.ne.0)then
                     24: c        write(0,1001) info,u,d1,d2,d3
                     25: c 1001   format(' sgefa abort. info=',i5,/4(3e12.4/))
                     26: c        stop 1
                     27:          if(d1(1).eq.d2(1)) then
                     28:             if(d1(2).eq.d2(2)) then
                     29:                u(3) = d1(3)
                     30:             else
                     31:                u(3) = d1(3)+(d2(3)-d1(3))*(u(2)-d1(2))/(d2(2)-d1(2))
                     32:             end if
                     33:          else
                     34:             u(3) = d1(3)+(d2(3)-d1(3))*(u(1)-d1(1))/(d2(1)-d1(1))
                     35:          end if
                     36:         return
                     37:       end if
                     38:       call sgesl(a,3,3,ipvt,r,0)
                     39:       u(3)=r(1)+r(2)*u(1)+r(3)*u(2)
                     40:       end

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.