home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume1 / 8712 / 2 < prev    next >
Encoding:
Internet Message Format  |  1990-07-13  |  11.1 KB

  1. Path: uunet!husc6!think!ames!sri-spam!rutgers!clyde!cbosgd!mandrill!hal!ncoast!allbery
  2. From: edward@pur-ee.UUCP (Edward L Haletky)
  3. Newsgroups: comp.sources.misc
  4. Subject: 2D contouring routines
  5. Message-ID: <6331@ncoast.UUCP>
  6. Date: 4 Dec 87 01:41:45 GMT
  7. Sender: allbery@ncoast.UUCP
  8. Organization: Purdue University Engineering Computer Network
  9. Lines: 437
  10. Approved: allbery@ncoast.UUCP
  11. X-Archive: comp.sources.misc/8712/2
  12.  
  13.  
  14.  
  15.  
  16. NetLanders;
  17.  
  18. The contour routine is not big. It is actually quite small.
  19. The key to using the routines is as follows:
  20.  
  21.         The routines expect the data to be upon some plane of coordinates x
  22. and  y or whatever is being used. The random points are located within this
  23. 2d mesh of x and y points. The routine requires that 3 2d arrays should  be
  24. present the first array is the location of the x point in mesh coordinates,
  25. and so on. You can read about it in the Man page supplied.  This  could  be
  26. used to do random point contouring by masking the unused coordinates within
  27. the arrays with some value that you know the random values will not  touch.
  28. Either way above the max or way below the minimun value of the function  to
  29. be  plotted.  The  source and man page follow. It uses unix type plot calls
  30. from a package developed here at Purdue (CRC  Graphics).  I  am  sure  your
  31. system has a similiar package....Ed 
  32. ---------------8<------------------------8<--------------------8<-------------
  33. # This is a shell archive.  Save this into a file, edit it
  34. # and delete all lines above this comment.  Then give this
  35. # file to sh by executing the command "sh file".  The files
  36. # will be extracted into the current directory owned by
  37. # you with default permissions.
  38. #
  39. # The files contained herein are:
  40. #      UContour.f    UContour.3G
  41. #
  42. echo 'x - UContour.f'
  43. sed 's/^X//' <<'________This_Is_The_END________' >>UContour.f
  44. Xc
  45. Xc    UContour was originally written by Dr. W.J. Usab for the
  46. Xc    PEPL graphics package. The routines have been modified to
  47. Xc    use the CRC graphics package by Edward L. Haletky and Tom 
  48. Xc    Jentink.
  49. Xc
  50. Xc    The following modifications were made: 11/13/87
  51. Xc        Use of output files and not direct printing (gf)
  52. Xc        gf options: filename - plot to file filename
  53. Xc                   scr      - plot to standard output
  54. Xc                   lpr      - plot to Printronix
  55. Xc
  56. Xc        added ml option - specifies length of column in 2D matrix
  57. Xc        Call to plots,fname were added to routines
  58. Xc        Added Label option
  59. Xc        Call to plot(0,0,999) added
  60. Xc        Number call changed to reflect C syntax on formats
  61. Xc        Transformed to lower case
  62. Xc
  63. Xc    To Compile:
  64. Xc        f77 -c UContour.f -lG
  65. Xc
  66. Xc
  67. Xc      subroutine contour(x,y,f,n,m,cmin,cmax,inc)
  68. X    subroutine ucontour(x,y,f,ml,in,jmax,cmin,cmax,inc,gf,
  69. X     +lbl,bd)
  70. X    real xmin,xmax,ymin,ymax,gscale
  71. X        common /climit/ xmin,xmax,ymin,ymax,gscale
  72. X    character *8 option
  73. X    character *20 lbl
  74. X    character *3 dev1,dev2,gf
  75. X    integer dev,blank
  76. X    real x(ml,*),y(ml,*),f(ml,*),bd(4),jmax(ml)
  77. Xc    define constants
  78. X        option = ' '
  79. X        dev1 = "scr"
  80. X        dev2 = "lpr"
  81. X    dev = 0
  82. X        if (gf.eq.dev1) dev = 6
  83. X        if (gf.eq.dev2) dev = 16
  84. X    blank = 0
  85. X    xfac  = 10.0
  86. X    yfac  = 10.0
  87. X
  88. Xc    define min and max ranges for x,y,f
  89. X
  90. X    call plots(dev,blank,option)
  91. X    if ((dev.ne.16).and.(dev.ne.6)) call fname(gf)
  92. X    call symbol(2,9.5,0.075,lbl,0.0)
  93. Xc
  94. Xc         define min and max ranges for x,y,f
  95. Xc
  96. X      xmin = x(1,1)
  97. X      xmax = x(1,1)
  98. X      ymin = y(1,1)
  99. X      ymax = y(1,1)
  100. X      fmin = f(1,1)
  101. X      fmax = f(1,1)
  102. Xc
  103. Xc      do 10 j = 1,m
  104. Xc      do 10 i = 1,n
  105. X      do 10 i=1,in
  106. X       m=jmax(i)
  107. X      do 10 j=1,m
  108. X      xmin = min(x(i,j),xmin)
  109. X      xmax = max(x(i,j),xmax)
  110. X      ymin = min(y(i,j),ymin)
  111. X      ymax = max(y(i,j),ymax)
  112. X      fmin = min(f(i,j),fmin)
  113. X      fmax = max(f(i,j),fmax)
  114. X   10 continue
  115. Xc
  116. X      xscale = 1.0*xfac/(xmax-xmin)
  117. X      yscale = 1.0*yfac/(ymax-ymin)
  118. X      gscale = min(xscale,yscale)
  119. Xc
  120. Xc      do 20 j = 1,m-1
  121. Xc      do 20 i = 1,n-1
  122. X      do 20 i=1,in-1
  123. X       m=jmax(i)
  124. X      do 20 j=1,m-1
  125. X      x1 = gscale*(x(i,j)-xmin)
  126. X      y1 = gscale*(y(i,j)-ymin)
  127. X      f1 = f(i,j)
  128. X      x2 = gscale*(x(i,j+1)-xmin)
  129. X      y2 = gscale*(y(i,j+1)-ymin)
  130. X      f2 = f(i,j+1)
  131. Xc******************************
  132. X       if (j.eq.(m-1)) then  
  133. X         x3=gscale*x(i+1,j)
  134. X         y3=gscale*y(i+1,j)
  135. X         f3=f(i+1,j)
  136. X        else
  137. X      x3 = gscale*(x(i+1,j+1)-xmin)
  138. X      y3 = gscale*(y(i+1,j+1)-ymin)
  139. X      f3 = f(i+1,j+1)
  140. X       endif
  141. Xc******************************
  142. X      x4 = gscale*(x(i+1,j)-xmin)
  143. X      y4 = gscale*(y(i+1,j)-ymin)
  144. X      f4 = f(i+1,j)
  145. X      diag1 = (x3-x1)**2+(y3-y1)**2
  146. X      diag2 = (x4-x2)**2+(y4-y2)**2
  147. X      if(diag1.le.diag2) then
  148. X        call tricon(x1,y1,f1,x2,y2,f2,x3,y3,f3,
  149. X     1              cmin,cmax,inc)
  150. X        call tricon(x1,y1,f1,x3,y3,f3,x4,y4,f4,
  151. X     1              cmin,cmax,inc)
  152. X      else
  153. X        call tricon(x1,y1,f1,x2,y2,f2,x4,y4,f4,
  154. X     1              cmin,cmax,inc)
  155. X        call tricon(x2,y2,f2,x3,y3,f3,x4,y4,f4,
  156. X     1              cmin,cmax,inc)
  157. X      endif
  158. Xc
  159. X   20 continue
  160. Xc
  161. Xc       plot contour region boundary (clockwise from 0,0)
  162. Xc
  163. X      x1 = gscale*(x(1,1)-xmin)
  164. X      y1 = gscale*(y(1,1)-ymin)
  165. X      call plot(x1,y1,3)
  166. Xc
  167. X       m=jmax(1)
  168. X       x1=gscale*x(1,m)
  169. X       y1=gscale*y(1,m)
  170. X       call plot(x1,y1,2)
  171. Xc
  172. X      
  173. X       x1=gscale*bd(1)      
  174. X       y1=gscale*bd(2)
  175. X       call  plot(x1,y1,2)
  176. Xc
  177. X        x1=gscale*bd(3)
  178. X        y1=gscale*bd(4)
  179. X        call  plot(x1,y1,2)
  180. Xc
  181. X       if (bd(3).lt.x(in,1)) then
  182. X       x1=gscale*x(in,1)  
  183. X       y1=gscale*bd(4)
  184. X       call  plot(x1,y1,2)
  185. Xc
  186. X       x1 = gscale*x(in,1)
  187. X       y1 = gscale*y(in,1)
  188. X       call plot(x1,y1,2)
  189. Xc
  190. X      x1 = gscale*(x(1,1)-xmin)
  191. X      y1 = gscale*(y(1,1)-ymin)
  192. X      call plot(x1,y1,2)
  193. X      endif
  194. X    call plot(0,0,999)
  195. Xc
  196. Xc
  197. X      return
  198. X      end
  199. Xc
  200. X      subroutine tricon(xa,ya,fa,xb,yb,fb,xc,yc,fc,cmin,cmax,inc)
  201. X      real x(3),y(3),f(3)
  202. Xc
  203. X      x(1) = xa
  204. X      x(2) = xb
  205. X      x(3) = xc
  206. X      y(1) = ya
  207. X      y(2) = yb
  208. X      y(3) = yc
  209. X      f(1) = fa
  210. X      f(2) = fb
  211. X      f(3) = fc
  212. Xc
  213. Xc         sort points in ascending order 
  214. Xc 
  215. X    do 10 n=1,2
  216. X    do 10 i=1,3-n
  217. X    if(f(i).gt.f(i+1))then
  218. X      xt = x(i) 
  219. X      yt = y(i)
  220. X      ft = f(i)
  221. X      x(i) = x(i+1)
  222. X      y(i) = y(i+1)
  223. X      f(i) = f(i+1)
  224. X      x(i+1) = xt
  225. X      y(i+1) = yt
  226. X      f(i+1) = ft
  227. X    endif
  228. X 10     continue
  229. Xc
  230. Xc         locate contour range within triangle
  231. Xc
  232. X      dc = (cmax-cmin)/float(inc)
  233. X      icmin = int((f(1)-cmin)/dc)
  234. X      icmin = max(icmin,0)
  235. X      icmax = int((f(3)-cmin)/dc)
  236. X      icmax = min(icmax,inc)
  237. X      if(icmax.lt.icmin) return
  238. Xc
  239. Xc         plot contours within range
  240. Xc
  241. X      do 20 i = icmin,icmax
  242. Xc       do 20 i=0,inc
  243. X       clev = cmin+dc*float(i)
  244. Xc      if(clev.lt.f(1)) go to 20
  245. X      if (clev.le.f(1)) goto 20
  246. Xc      if(clev.gt.f(3)) go to 20
  247. X      if (clev.ge.f(3)) goto 20
  248. X      if(f(2).gt.clev) then
  249. X        fac  = (clev-f(1))/(f(2)-f(1))
  250. X        x1   = x(1)+fac*(x(2)-x(1))
  251. X        y1   = y(1)+fac*(y(2)-y(1))
  252. X      else
  253. X        fac  = (clev-f(2))/(f(3)-f(2))
  254. X        x1   = x(2)+fac*(x(3)-x(2))
  255. X        y1   = y(2)+fac*(y(3)-y(2))
  256. X      endif
  257. X      fac  = (clev-f(1))/(f(3)-f(1))
  258. X      x2   = x(1)+fac*(x(3)-x(1))
  259. X      y2   = y(1)+fac*(y(3)-y(1))
  260. Xc
  261. X      call plot(x1,y1,3)
  262. X      call plot(x2,y2,2)
  263. Xc      call labinc(i,clev,x1,y1,x2,y2)
  264. X   20 continue
  265. Xc
  266. X      return
  267. X      end
  268. Xc
  269. X      subroutine labinc(i,clev,x1,y1,x2,y2)
  270. Xc      common /clabinc/ icount
  271. Xc      data icount/0/
  272. X      tdist = abs(x2-x1)+abs(y2-y1)
  273. X      if(tdist.lt.1e-20) then
  274. X        write(6,*)' ** error in labinc: tdist=',tdist
  275. X        return
  276. X      endif
  277. X      interv = 21
  278. X      icount = icount+1
  279. X      if(icount.lt.interv) return
  280. X      icount = 0
  281. X      pi = acos(-1.0)
  282. X      hpi = 0.5*pi
  283. X      theta = atan2((y2-y1),(x2-x1))
  284. X      if(theta.ge.hpi) then
  285. X        theta = theta-pi
  286. X      else if(theta.lt.-hpi) then
  287. X        theta = theta+pi
  288. X      endif
  289. X      thetan = theta+hpi
  290. X      x   = 0.5*(x1+x2)+0.05*cos(thetan)-0.225*cos(theta)
  291. X      y   = 0.5*(y1+y2)+0.05*sin(thetan)-0.225*sin(theta)
  292. X      thetad = 180.0*theta/pi
  293. Xcc      call number(x,y,.1,i,0.0,'i2')
  294. Xc      call number(x,y,.075,clev,thetad,'f6.3')
  295. Xc       call number(x,y,.055,thetad,"%6.3f",clev)
  296. X      return
  297. X      end
  298. ________This_Is_The_END________
  299. echo 'x - UContour.3G'
  300. sed 's/^X//' <<'________This_Is_The_END________' >>UContour.3G
  301. X.\" troff -man
  302. X.EN
  303. X.TH UContour 3G  11/13/87
  304. X.SH NAME
  305. XUContour - EGraphics Contour Routine.
  306. X.SH AUTHORS
  307. XDr. W.J.Usab Jr.
  308. X.br
  309. XEdward L. Haletky
  310. X.br
  311. XTom Jentink
  312. X.SH SYNOPSIS
  313. Xcall ucontour(x,y,f,ml,in,jmax,cmin,cmax,inc,gf,lbl,bd)
  314. X.SH DESCRIPTION
  315. X.fi
  316. XUContour will plot a family of curves of number 
  317. X.I inc
  318. Xwithin the range
  319. X.I [cmin:cmax].
  320. XAll the curves will be of a constant value of the function,
  321. X.I f.
  322. X.sp
  323. X.RS
  324. X.fi
  325. X.B x
  326. X- two dimensional array for x-coordinate at grid point (n,m) at which
  327. Xthe function, 
  328. X.I f,
  329. Xwas computed.
  330. X.sp 1
  331. X.B y
  332. X- two dimensional array for y-coordinate at grid point (n,m) at which
  333. Xthe function, 
  334. X.I f,
  335. Xwas computed.
  336. X.sp 1
  337. X.B f 
  338. X- two dimensional array for the function at grid point (n,m).
  339. X.sp 1
  340. X.B ml
  341. X- maximun length of the arrays
  342. X.I x,
  343. X.I y,
  344. X.I f.
  345. XThe arrays are declared in the following manor:
  346. X.br
  347. X.ce 1
  348. Xreal x(ml,*),y(ml,*),f(ml,*)
  349. X.br
  350. X.sp 1
  351. X.B in 
  352. X- maximun number of grid points in n-direction.
  353. X.sp 1
  354. X.B jmax  
  355. X- an array that gives maximun number of grid points in the m-direction,
  356. Xfor each,
  357. X.I n
  358. Xstation.
  359. X.sp 1
  360. X.B cmin 
  361. X- minimun value of the constant-value contour of the function,
  362. X.I f.
  363. X.sp 1
  364. X.B cmax 
  365. X- maximun value of the constant-value contour of the function,
  366. X.I f.
  367. X.sp 1
  368. X.B inc 
  369. X- integer number of contours between the range
  370. X.I [cmin:cmax].
  371. X.sp 1
  372. X.B gf  
  373. X- device which to plot
  374. X.RE
  375. X.in +10
  376. X.B scr
  377. X- plot to stdout.
  378. X.br
  379. X.B lpr
  380. X- plot to Printronix Printer.
  381. X.br
  382. X.B gf
  383. X- plot to file. File can be printed with the
  384. X.I gplp(1G)
  385. Xcommand. The 
  386. X.I filename
  387. Xcan not be more than 3 characters long.
  388. X.in -10
  389. X.RS
  390. X.sp 1
  391. X.B lbl
  392. X- label the contour plot.
  393. X.B lbl
  394. Xcan be a maximun of 20 characters long.
  395. X.br
  396. X.sp 1
  397. X.B bd
  398. X- an array that holds the bounding area values.
  399. X.RE
  400. X.in 10
  401. X.B bd(1)
  402. X- xb1
  403. X.br
  404. X.B bd(2)
  405. X- yb1
  406. X.br
  407. X.B bd(3)
  408. X- xb2
  409. X.br
  410. X.B bd(4)
  411. X- yb2
  412. X.br
  413. X.in -10
  414. X.RS
  415. X.sp 1
  416. X.B "To Compile"
  417. Xf77 prog.f -L/x/edward/lib -lE -lG
  418. X.RE
  419. X.SH FILES
  420. X.RS
  421. X/x/edward/lib/libE.a
  422. X.br
  423. X/x/edward/src/UContour.f
  424. X.RE
  425. X.SH DIAGNOSTICS
  426. X.LP
  427. XUContour will blank the page before plotting. No CRC graphics options are used.
  428. XThe Border of the Channel will be printed, it is a crude representation. The
  429. Xlabeling of the lines will be implemented at a later date. There are a few bugs
  430. Xstill.
  431. X.SH "SEE ALSO"
  432. X.B "CRC Graphics Manaul ,"
  433. X.B gplp(1G)
  434. X.SH BUGS
  435. X.RS
  436. XPlease Report all bugs to:
  437. X.br
  438. Xedward@ga.ecn.purdue.edu
  439. X.RE
  440. ________This_Is_The_END________
  441. exit
  442. ---------------8<------------------------8<--------------------8<-------------
  443.  
  444.  
  445. =====================================    ^    =================================
  446.     Edward L. Haletky            |E     |o|     U| "To race against the sky..."
  447. Usenet: ~!pur-ee!edward          |L   ^/| |\^   S| "The world of the
  448. Arpa:   edward@ga.ecn.purdue.edu |H   / | | \   A|  immagination is boundless."
  449. =====================================/__-|-__\=================================
  450.