![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to box.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / view2d / Tri|
Return to box.f CVS log | Up to [Research Unix] / researchv10no / cmd / view2d / Tri |
1.1 root 1: ******************from mhuxt!rksmith on 25 Feb 85******************
2: *DECK INTBOX
3: subroutine intbox (nx,xmin,xmax,ny,ymin,ymax,f,tr,
4: * nv,fk,floor,vx,vy,itri)
5: dimension f(nx,1),tr(nx,1),fk(1),vx(1),vy(1),itri(3,1)
6: c
7: c intbox interpolates points from triangular grid
8: c (fk,vx,vy) onto rectangular grid ( f(nx,ny) )
9: c
10: c tr triangle indices
11: c f values on grid
12: c floor out-of-bounds
13: c itri(3,nt)
14: if (nx.lt.1.or.ny.lt.1) return
15: if (nx.eq.1) then
16: dx=1.0
17: else
18: dx=(xmax-xmin)/float(nx-1)
19: end if
20: c
21: if (ny.eq.1) then
22: dy=1.0
23: else
24: dy=(ymax-ymin)/float(ny-1)
25: end if
26: c
27: c The main loop
28: c
29: yp=ymin
30: do 200 iy=1,ny
31: xp=xmin
32: do 100 ix=1,nx
33: it=int(tr(ix,iy))
34: if (it.lt.1) then
35: f(ix,iy)=floor
36: else
37: c
38: c linear interpolation on triangle it
39: c
40: iv1=itri(1,it)
41: iv2=itri(2,it)
42: iv3=itri(3,it)
43: c
44: x13=vx(iv1)-vx(iv3)
45: x23=vx(iv2)-vx(iv3)
46: y13=vy(iv1)-vy(iv3)
47: y23=vy(iv2)-vy(iv3)
48: det=x13*y23-y13*x23
49: c1=y23*(xp-vx(iv3))-x23*(yp-vy(iv3))
50: c2=x13*(yp-vy(iv3))-y13*(xp-vx(iv3))
51: c
52: df=c1*(fk(iv1)-fk(iv3))+c2*(fk(iv2)-fk(iv3))
53: f(ix,iy)=fk(iv3)+df/det
54: c
55: end if
56: 100 xp=xp+dx
57: 200 yp=yp+dy
58: c
59: return
60: end
61: *DECK MACBOX
62: subroutine macbox (nv,vx,vy,nt,itri,
63: * nx,xmin,xmax,ny,ymin,ymax,tr)
64: dimension vx(1),vy(1),itri(3,1)
65: dimension tr(nx,1)
66: c
67: c find triangles containing box points
68: c
69: if (nx.lt.1.or.ny.lt.1) return
70: if (nx.eq.1) then
71: dx=1.0
72: else
73: dx=(xmax-xmin)/float(nx-1)
74: if (dx.eq.0.0) dx=1.0
75: end if
76: dx1=1.0/dx
77: xdx=xmin-dx
78: c
79: if (ny.eq.1) then
80: dy=1.0
81: else
82: dy=(ymax-ymin)/float(ny-1)
83: if (dy.eq.0.0) dy=1.0
84: end if
85: dy1=1.0/dy
86: ydy=ymin-dy
87: c
88: c search triangles for points
89: c
90: eps=1.0e-5
91: cnorm=1.0-4.0*eps
92: do 202 iy=1,ny
93: do 201 ix=1,nx
94: tr(ix,iy)=0.0
95: 201 continue
96: 202 continue
97: c
98: do 230 it=1,nt
99: c
100: iv1=itri(1,it)
101: iv2=itri(2,it)
102: iv3=itri(3,it)
103: c
104: vxmin=dx1*(amin1(vx(iv1),vx(iv2),vx(iv3))-xmin)+1.0
105: ixmin=max0(int(vxmin),1)
106: vxmax=dx1*(amax1(vx(iv1),vx(iv2),vx(iv3))-xmin)+1.0
107: ixmax=min0(int(vxmax),nx)
108: vymin=dy1*(amin1(vy(iv1),vy(iv2),vy(iv3))-ymin)+1.0
109: iymin=max0(int(vymin),1)
110: vymax=dy1*(amax1(vy(iv1),vy(iv2),vy(iv3))-ymin)+1.0
111: iymax=min0(int(vymax),ny)
112: c
113: c find box points contained in this macro triangle
114: c
115: x13=vx(iv1)-vx(iv3)
116: x23=vx(iv2)-vx(iv3)
117: y13=vy(iv1)-vy(iv3)
118: y23=vy(iv2)-vy(iv3)
119: det=x13*y23-y13*x23
120: c
121: x3=xdx-vx(iv3)
122: y3=ydy-vy(iv3)
123: c
124: do 220 iy=iymin,iymax
125: yp=y3+float(iy)*dy
126: do 210 ix=ixmin,ixmax
127: xp=x3+float(ix)*dx
128: c
129: c1=(y23*xp-x23*yp)/det
130: if (c1+eps.lt.0.0) go to 210
131: c2=(x13*yp-y13*xp)/det
132: if (c2+eps.lt.0.0) go to 210
133: c3=1.0-c1-c2
134: if (c3+eps.lt.0.0) go to 210
135: c
136: tr(ix,iy)=float(it)+0.5
137: c
138: 210 continue
139: 220 continue
140: c
141: 230 continue
142: c
143: return
144: end
145:
146:
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.