00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 subroutine plt3d(a,m,n,alt,az,xlen,xoff,ylen,yoff,zfac,zoff,ier)
00021 save
00022 c
00023 c--- 3d plot routine, adapted from the SAS-3 routine of the same name.
00024 c
00025 c--- calling sequence:
00026 c call plt3d (a,m,n,alt,az,xlen,xoff,ylen,yoff,zfac,zoff,ier)
00027 c
00028 c a - real data array, represents height of surface as
00029 c function of location in plane
00030 c m,n - dimensions of data array a
00031 c alt, az - altitude, azimuth viewing angles in degrees
00032 c xlen,ylen - length of unprojected x,y axes in inches
00033 c xoff,yoff - offset of origin in inches
00034 c zfac - scaling of z-axis from data units to unprojected inches
00035 c zoff - offset of z-origin in user units
00036 c ier - returns 0 for no error
00037 c
00038 c Altitude is measured from the x-y plane toward the positive z-axis.
00039 c Azimuth is measured around the z-axis (as seen from above),
00040 c counterclockwise from the positive x-axis.
00041 c
00042 c The origin (xoff, yoff), in screen coordinates, is the image of the
00043 c point at the CENTER of the x-y grid and at a height zoff above it
00044 c (in the units used for a).
00045 c
00046 dimension a(m,n)
00047 c
00048 c Note that the original user-provided work arrays are now internal
00049 c and static.
00050 c
00051 parameter (LDIM = 5000)
00052 dimension x(LDIM),y(LDIM)
00053 c
00054 common /plt3b/ a1,a2,a3,b1,b2,b3,b4,zoff1
00055 c
00056 parameter (PI = 3.141592654, PIDEG = PI/180.)
00057 c
00058 c Coordinate to plot conversions:
00059 c
00060 c (1,1) <--> (xmin,ymin)
00061 c (1,n) <--> (xmin,ymax)
00062 c (m,1) <--> (xmax,ymin)
00063 c (m,n) <--> (xmax,ymax)
00064 c
00065 xplot(i,j) = a1*i + a2*j + a3
00066 yplot(i,j) = b1*i + b2*j + b3*(a(i,j)-zoff1) + b4
00067 c
00068 lmax = 2*min(m,n)
00069 if (LDIM.lt.lmax) then
00070 write(0,*)'PLT3D: Array too big!'
00071 ier = 2
00072 return
00073 end if
00074 c
00075 taz = az*PIDEG
00076 talt = alt*PIDEG
00077 c
00078 saz = sin(taz)
00079 caz = cos(taz)
00080 sal = sin(talt)
00081 cal = cos(talt)
00082 c
00083 xsc = xlen/float(m-1)
00084 ysc = ylen/float(n-1)
00085 c
00086 c The viewing plane is such that the "horizontal" vector is
00087 c
00088 c x' = (-saz, caz, 0),
00089 c
00090 c and the "vertical" vector is
00091 c
00092 c y' = (-sal*caz, -sal*saz, cal).
00093 c
00094 a1 = -saz*xsc
00095 a2 = caz*ysc
00096 a3 = xoff - 0.5*(a1*float(m+1)+a2*float(n+1))
00097 b1 = -caz*sal*xsc
00098 b2 = -saz*sal*ysc
00099 b3 = zfac*cal
00100 b4 = yoff - 0.5*(b1*float(m+1)+b2*float(n+1))
00101 zoff1 = zoff
00102 c
00103 c Find which quadrant we are viewing from.
00104 c
00105 iaz = 1
00106 if (caz.le.0.0) iaz = iaz + 1
00107 if (saz.le.0.0) iaz = 5 - iaz
00108 c
00109 c Choose the direction of traversal of the plot based on view angle.
00110 c
00111 if (iaz.eq.2.or.iaz.eq.3) then
00112 ifirst = 1
00113 else
00114 ifirst = m
00115 end if
00116 ilast = m + 1 - ifirst
00117 istep = sign(1,ilast-ifirst)
00118 c
00119 if (iaz.ge.3) then
00120 jfirst = 1
00121 else
00122 jfirst = n
00123 end if
00124 jlast = n + 1 - jfirst
00125 jstep = sign(1,jlast-jfirst)
00126 c
00127 c Make sure the array sent to nxtvu is traversed from left to right.
00128 c
00129 if (iaz.eq.1.or.iaz.eq.3) then
00130 lli = -1
00131 else
00132 lli = 1
00133 end if
00134 c
00135 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00136 c
00137 ic = 0
00138 ibeg = ifirst+istep
00139 c
00140 c Traverse the plot region in a zig-zag motion, transverse to the
00141 c viewing direction (left to right), and from front to back.
00142 c The plotting routine "nxtvu" takes care of hidden line removal.
00143 c
00144 c First, deal with the "front" portion of the plot (i.e. nearest
00145 c the viewer).
00146 c
00147 70 lnth = min(2*iabs(ibeg-ifirst)+1,lmax)
00148 if (lli.eq.-1) then
00149 ll = lnth + 1
00150 else
00151 ll = 0
00152 end if
00153 c
00154 i = ibeg
00155 j = jfirst
00156 ll = ll+lli
00157 x(ll) = xplot(i,j)
00158 y(ll) = yplot(i,j)
00159 c
00160 80 i = i-istep
00161 ll = ll+lli
00162 x(ll) = xplot(i,j)
00163 y(ll) = yplot(i,j)
00164 if (j.ne.jlast) then
00165 j = j+jstep
00166 ll = ll+lli
00167 x(ll) = xplot(i,j)
00168 y(ll) = yplot(i,j)
00169 if (i.ne.ifirst) go to 80
00170 end if
00171 c
00172 call nxtvu(ic,x,y,lnth,ier)
00173 if (ier.ne.0) return
00174 ic = 1
00175 if (ibeg.ne.ilast) then
00176 ibeg = ibeg+istep
00177 go to 70
00178 end if
00179 c
00180 c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
00181 c
00182 c Now do the "rear" half of the figure in the same manner.
00183 c
00184 jbeg = jfirst
00185 c
00186 100 lnth = min(2*iabs(jbeg-jlast)+1,lmax)
00187 if (lli.eq.-1) then
00188 ll = lnth + 1
00189 else
00190 ll = 0
00191 end if
00192 c
00193 i = ilast
00194 j = jbeg
00195 ll = ll+lli
00196 x(ll) = xplot(i,j)
00197 y(ll) = yplot(i,j)
00198 c
00199 110 j = j+jstep
00200 ll = ll+lli
00201 x(ll) = xplot(i,j)
00202 y(ll) = yplot(i,j)
00203 if (i.ne.ifirst) then
00204 i = i-istep
00205 ll = ll+lli
00206 x(ll) = xplot(i,j)
00207 y(ll) = yplot(i,j)
00208 if (j.ne.jlast) go to 110
00209 end if
00210 c
00211 call nxtvu(1,x,y,lnth,ier)
00212 if (ier.ne.0) return
00213 jbeg = jbeg+jstep
00214 if (jbeg.eq.jlast) return
00215 go to 100
00216 c
00217 end
00218
00219
00220 function xplot3d(x,y,z)
00221 save
00222 c
00223 c Externally accessible version of the statement function in plt3d.
00224 c This will work only AFTER plt3d has been used, as that routine
00225 c sets up the common array.
00226 c
00227 c The coordinates x and y are expected to range from 1 to m and
00228 c 1 to n, respectively.
00229 c
00230 c Typical usage: xx = xplot3d(float(i),float(j),array(i,j))
00231 c
00232 c Note that z is not used.
00233 c
00234 common /plt3b/ a1,a2,a3,b1,b2,b3,b4,zoff1
00235 c
00236 xplot3d = a1*x + a2*y + a3
00237 c
00238 end
00239
00240
00241 function yplot3d(x,y,z)
00242 save
00243 c
00244 c Externally accessible version of the statement function in plt3d.
00245 c This will work only AFTER plt3d has been used, as that routine
00246 c sets up the common array.
00247 c
00248 c The coordinates x and y are expected to range from 1 to m and
00249 c 1 to n, respectively.
00250 c
00251 common /plt3b/ a1,a2,a3,b1,b2,b3,b4,zoff1
00252 c
00253 c Typical usage: yy = yplot3d(float(i),float(j),array(i,j))
00254 c
00255 yplot3d = b1*x + b2*y + b3*(z-zoff1) + b4
00256 c
00257 end
