|
|
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.