# farb2d - Fill area with bicubics for 2d contour plotting # # Description # This routine produces contours of the given data on the graphics # output. # # This code is based on the ACM algorithm 671 published in the # _ACM_Transactions_of_Mathematical_Software, Vol 15, No. 1, March 1989, # Pages 79-89. The source code itself was retrieved by the Internet mail # service netlib located at ornl.gov. For more information on netlib, send # mail with a single line containing "help" to netlib@ornl.gov. # # NOTE: BESIDES THE COPYRIGHTS OF STSCI/STSDAS, THIS SOFTWARE IS SUBJECT # TO COPYRIGHTS OF THE ACM. CONTACT THE ACM FOR MORE INFORMATION. # # Below this introduction is the header appended by netlib. Below that # and for each subroutine, the header by the author is also included. # Besides the rewrite into SPP and the inclusion of the IRAF GIO graphics # interface, the code has not been modified. This includes the preservation # of GOTO's. Only in the simplest cases have GOTO's been replaced with # other constructs. # # NOTE: A logical change has been made. The problem is with the routine # usrplt, specifically the use of the ncol parameter. The documentation # seemed very clear on the point- ncol is the color (or color index into # some color table) to draw the area or line in. However, in usage, ncol # seems to have two distinct usages. If the mode to draw is fill area, # the usage is consistent with documentation. However, if the mode is # line drawing only, ncol was not set to the color, but was given the # contour level index (the value representing which contour level was # being drawn). This is so highly inconsistent that it gives one a # massive head trauma thinking about it. The solution here is to change # the calls involving line drawing mode to pass the color index and not # the contour index. # # History # 31Dec90 - Retrieved from netlib. Prototype for applicability to IRAF. # Jonathan D. Eisenhamer, STScI. # 15Jan91 - Rewritten into SPP. jde # 21Jan91 - Modified selected calls to usrplt to make the usrplt interface # consistent. jde # # Bugs # - Though these routines handle both simple contouring by lines and # area filling, only the line drawing works. The IRAF GIO gfill # area filling is not implemented (insert appropriate comments by the # reader). However, if gfill is ever implemented, nothing in these # routines should need modifying except to debug the graphics calls # for correct usage. # #--------------------------------------------------------------------------- # Intro added by netlib: #--------------------------------------------------------------------------- # From netlibd@surfer.EPM.ORNL.GOV Mon Dec 31 08:19:16 1990 # To: eisenham@stsci.edu # Subject: send 671 from toms # # Disclaimer. This on-line distribution of algorithms is not an official # service of the Association for Computing Machinery (ACM). This service is # being provided on a trial basis and may be discontinued in the future. ACM # is not responsible for any transmission errors in algorithms received by # on-line distribution. Official ACM versions of algorithms may be obtained # from: IMSL, Inc. # 2500 ParkWest Tower One # 2500 CityWest Blvd. # Houston, TX 77042-3020 # (713) 782-6060 # # Columns 73-80 have been deleted, trailing blanks removed, and program # text mapped to lower case. Assembly language programs have sometimes # been deleted and standard machine constant and error handling # routines removed from the individual algorithms and collected # into the libraries "core" and "slatec". In a few cases, this trimming # was overzealous, and will be repaired when we can get back the original # files. We have tried to incorporate published Remarks; if we missed # something, please contact Eric Grosse, 201-582-5828, research!ehg or # ehg@research.att.com. # # The material in this library is copyrighted by the ACM, which grants # general permission to distribute provided the copies are not made for # direct commercial advantage. For details of the copyright and # dissemination agreement, consult the current issue of TOMS. # # Caveat receptor. (Jack) dongarra@cs.utk.edu, (Eric Grosse) research!ehg # Careful! Anything free comes with no guarantee. # * # *We would like to thank Sequent Computer Systems Inc. # *for the computer used to run netlib. # * # *** from netlib, Mon Dec 31 08:20:42 EST 1990 *** # C ALGORITHM 671, COLLECTED ALGORITHMS FROM ACM. # C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, # C VOL. 15, NO. 1, PP. 79-89. # FARB-E-2D Version 2.1 10/1988 # #--------------------------------------------------------------------------- # End of netlib intro. # Author intro: #--------------------------------------------------------------------------- # # fill area with bicubics for 2d contour plotting # ----------------------------------------------- # farb-e-2d version 2.1, 10/1988 # # t r i p algorithm a. preusser # # author: a. preusser # fritz-haber-institut der mpg # faradayweg 4-6 # d-1000 berlin 33 # # input parameters # x array of length lx for x-coordinates of # a regular grid # in ascending order. # # x- and y-coordinates must be given # in centimeters # ============== # # lx number of grid lines x= x(i), i=1,lx # parallel to y-axis. # y array of length ly for y-coordinates # in ascending order. # ly number of grid lines y= y(i), i=1,ly # parallel to x-axis. # z 2-dimensional array dimensioned z(nxdim,...) # defining the z-values at the grid points. # the point with the coordinates x(k), y(l) # receives the value z(k,l), k=1,lx, l=1,ly. # nxdim first dimension of array z # cn array of length nc for the z-values of # the contours (contour levels) # in ascending order # icol integer array of length nc+1 for # the colours to be used for the lines or areas. # values from this array are passed to # the user supplied subroutine usrplt. # icol(i) is used for the area, where # z > cn(i-1) and # z <= cn(i), # for i=2,nc. # areas, where z <= cn(1) # are filled with colour icol(1), # and areas, where z > icol(nc) # are filled with colour icol(nc+1). # nc number of contour levels, nc <= 100 # mode 0, fill area only # 1, lines only # 2, fill area and lines # # # output # is performed by calls to the subroutine usrplt # to be supplied by the user (an example for usrplt # is included.) # # parameters of usrplt # subroutine usrplt (x,y,n,ncol,mode) # x,y real arrays of length n for # the coordinates of the polygon # to be plotted. # n number of points of polygon # ncol colour to be used # for the area or the line. # for ncol, the program passes # values of icol as described above. # mode 1, line drawing # 0, fill area # ------------------------------------------------------------- # # this module (farb2d) is based on subroutine sfcfit of # acm algorithm 474 by h.akima #--------------------------------------------------------------------------- # End of author intro. Now to code. #--------------------------------------------------------------------------- # Define the generic error code. define FATAL_ERROR 1 # Define some labels for gotos. define label_15_ 15 define label_16_ 16 define label_17_ 17 define label_18_ 18 define label_19_ 19 define label_20_ 20 define label_21_ 21 define label_25_ 25 define label_26_ 26 define label_27_ 27 define label_28_ 28 define label_36_ 36 define label_39_ 39 define label_40_ 40 define label_49_ 49 define label_56_ 56 define label_67_ 67 define label_68_ 68 define label_69_ 69 define label_85_ 85 define label_90_ 90 define label_30_ 30 procedure farb2d( x,lx,y,ly,z,nxdim,cn,icol,nc,mode ) real x[lx] # Array of length lx for x-coordinates of # a regular grid in ascending order. int lx # Length of array x. real y[ly] # Array of length ly for y-coordinates of # a regular grid in ascending order. int ly # Length of array y. real z[nxdim, ly] # 2-dimensional array dimensioned z(nxdim,...] # defining the z-values at the grid points. # The point with the coordinates x(k], y(l] # receives the value z(k,l], k=1,lx, l=1,ly. int nxdim # First dimension of array z. real cn[nc] # Array of length nc for the z-values of # the contours (contour levels) in ascending order. int icol[nc+1] # Integer array of length nc+1 for # the colours to be used for the lines or areas. # values from this array are passed to # the user supplied subroutine usrplt. # icol(i] is used for the area, where # z > cn(i-1] and # z <= cn(i], # for i=2,nc. areas, where z <= cn(1] # are filled with colour icol(1], # and areas, where z > icol(nc] # are filled with colour icol(nc+1]. int nc # Number of contour levels, nc <= 100. int mode # The mode to produce the plots in. Possible values: # 0 -> fill area only # 1 -> lines only # 2 -> fill area and lines # Declarations. # [jde: Sorry about no comments, but there weren't any in the original] real a5 real b1, b2, b3, b4, b5 real sw real w2, w3, wx2, wx3, wy2, wy3 real x3, x4, x5, x6, xx[4] real y3, y4, yy[4] real z33, z43, z53, z54, z5b1, z5b2, z5b3, z5b4, z5b5, z62, z65, z63, z64, z6b1, z6b2, z6b3, z6b4, z6b5, za[4,2], za5b2, za5b3, za5b4, zab[2,3], zb[5], zx[2], zx33, zx34, zxy[2], zxy33, zxy34, zy[2], zy33, zy34, zz[4], zzx[4], zzxy[4], zzy[4] int i, ifa, ix, ix6, ixm1, ixml, iy, iym2, iym3, iyml, iyml1 int jy int nsides int lx0, lxm1, lxm2, ly0, lym1, lym2 data x4,z43,x5,z53,z63,z5b1,z5b2,z5b3,z5b4,z5b5, a5,za5b2,za5b3,za5b4,x6,z64,z6b1,z6b2,z6b3,z6b4,z6b5, zx,zy,zxy,zab,za,zb /46*0/ begin # Preliminary processing. ifa= 0 lx0= lx lxm1= lx0-1 lxm2= lxm1 - 1 ly0= ly lym1= ly0 -1 lym2= lym1-1 # Error check if( lxm2 < 0 ) call error( FATAL_ERROR, "farb2d: x array is of 0 length" ) if( lym2 < 0 ) call error( FATAL_ERROR, "farb2d: y array is of 0 length" ) for ( ix = 2; ix <= lx0; ix = ix + 1 ) { if ( x[ix-1] - x[ix] == 0. ) call error( FATAL_ERROR, "farb2d: identical x values" ) if ( x[ix-1] - x[ix] > 0. ) call error( FATAL_ERROR, "farb2d: x values out of sequence" ) } for ( iy = 2; iy <= ly0; iy = iy + 1 ) { if ( y[iy-1] - y[iy] == 0. ) call error( FATAL_ERROR, "farb2d: identical y values" ) if ( y[iy-1] - y[iy] > 0. ) call error( FATAL_ERROR, "farb2d: y values out of sequence" ) } for ( i = 2; i <= nc; i = i + 1 ) { if ( cn[i-1] - cn[i] == 0. ) call error( FATAL_ERROR, "farb2d: identical contour values" ) if ( cn[i-1] - cn[i] > 0. ) call error( FATAL_ERROR, "farb2d: contour values out of sequence" ) } if( nc <= 0 ) call error( FATAL_ERROR, "farb2d: insufficient number of contours" ) # Main do-loops do iy=2,ly0 { iym2= iy-2 iym3= iy-3 iyml= iy-ly0 iyml1= iyml+1 ix6= 0 do ix=1,lx0 { ixm1= ix-1 ixml= ix- lx0 # routines to pick up necessary x,y, and z values to # compute the za,zb, and zab values, and to estimate # them when necessary # preliminary when ix == 1 if (ixm1 != 0) goto label_15_ y3= y[iy-1] y4= y[iy] b3= 1./(y4-y3) if (iym2 > 0) b2= 1./(y3-y[iy-2]) if (iym3 > 0) b1= 1./(y[iy-2]-y[iy-3]) if (iyml < 0) b4= 1./(y[iy+1]-y4) if (iyml1 < 0) b5= 1./(y[iy+2]-y[iy+1]) goto label_18_ # # to save the old values # label_15_ za[1,1]= za[2,1] za[1,2]= za[2,2] x3= x4 z33= z43 za[2,1]= za[3,1] za[2,2]= za[3,2] zab[1,1]= zab[2,1] zab[1,2]= zab[2,2] zab[1,3]= zab[2,3] label_16_ x4= x5 z43= z53 zb[1]= z5b1 zb[2]= z5b2 zb[3]= z5b3 zb[4]= z5b4 zb[5]= z5b5 za[3,1]= za[4,1] za[3,2]= za[4,2] zab[2,1]= za5b2 zab[2,2]= za5b3 zab[2,3]= za5b4 label_17_ x5= x6 z53= z63 z54= z64 z5b1= z6b1 z5b2= z6b2 z5b3= z6b3 z5b4= z6b4 z5b5= z6b5 # to compute the za, zb, and zab values and # to estimate the zb values # when (iy <= 3).or.(iy >= ly-1) label_18_ ix6= ix6 + 1 if (ix6 > lx0) goto label_26_ x6= x[ix6] z63= z[ix6,iy-1] z64= z[ix6,iy] z6b3= (z64-z63)*b3 if (lym2 == 0) goto label_20_ if (iym2 == 0) goto label_19_ z62= z[ix6,iy-2] z6b2= (z63-z62)*b2 if (iyml != 0) goto label_19_ z6b4= z6b3 + z6b3 -z6b2 goto label_21_ label_19_ z65= z[ix6,iy+1] z6b4= (z65-z64)*b4 if (iym2 != 0) goto label_21_ z6b2= z6b3 + z6b3 -z6b4 goto label_21_ label_20_ z6b2= z6b3 z6b4= z6b3 label_21_ if (iym3 > 0) z6b1= (z62-z[ix6,iy-3])*b1 else z6b1= z6b2 + z6b2 -z6b3 if (iyml1 < 0) z6b5= (z[ix6,iy+2] - z65)*b5 else z6b5= z6b4 + z6b4 -z6b3 if (ix6 == 1) goto label_17_ a5= 1./(x6-x5) za[4,1]= (z63-z53)*a5 za[4,2]= (z64-z54)*a5 za5b2= (z6b2-z5b2)*a5 za5b3= (z6b3-z5b3)*a5 za5b4= (z6b4-z5b4)*a5 if (ix6 == 2) goto label_16_ goto label_28_ # To estimate the za and zab values when (ix >= lx-1).and.(lx > 2) label_26_ if (lxm2 == 0) goto label_27_ za[4,1]= za[3,1] +za[3,1] - za[2,1] za[4,2]= za[3,2] +za[3,2] - za[2,2] if (ixml == 0) goto label_90_ za5b2= zab[2,1] + zab[2,1] - zab[1,1] za5b3= zab[2,2] + zab[2,2] - zab[1,2] za5b4= zab[2,3] + zab[2,3] - zab[1,3] goto label_90_ # To estimate the za and zab values when (ix >= lx-1).and.(lx == 2) label_27_ za[4,1]= za[3,1] za[4,2]= za[3,2] if (ixml == 0) goto label_90_ za5b2= zab[2,1] za5b3= zab[2,2] za5b4= zab[2,3] # To estimate the za and zab values when ix eq 1 label_28_ if (ixm1 != 0) goto label_90_ za[2,1]= za[3,1] +za[3,1] - za[4,1] za[1,1]= za[2,1] +za[2,1] - za[3,1] za[2,2]= za[3,2] +za[3,2] - za[4,2] za[1,2]= za[2,2] +za[2,2] - za[3,2] zab[1,1]= zab[2,1] + zab[2,1] -za5b2 zab[1,2]= zab[2,2] + zab[2,2] -za5b3 zab[1,3]= zab[2,3] + zab[2,3] -za5b4 goto label_30_ # Numerical differentation ---- to determine # partial deriv. zx,zy, and zxy as weighted means of # divided differences za, zb, and zab, respectively # to save the old values when ix != 1 label_90_ zx33= zx[1] zx34= zx[2] zy33= zy[1] zy34= zy[2] zxy33= zxy[1] zxy34= zxy[2] # # new computation label_30_ do jy=1,2 { w2= abs(za[4,jy]-za[3,jy]) w3= abs(za[2,jy]-za[1,jy]) sw= w2 + w3 if (sw != 0) { wx2= w2/sw wx3= w3/sw } else { wx2= 0.5 wx3= 0.5 } zx[jy]= wx2*za[2,jy] + wx3*za[3,jy] w2= abs(zb[jy+3]-zb[jy+2]) w3= abs(zb[jy+1]-zb[jy]) sw= w2 + w3 if (sw != 0.) { wy2= w2/sw wy3= w3/sw } else { wy2= 0.5 wy3= 0.5 } zy[jy]= wy2*zb[jy+1] + wy3*zb[jy+2] zxy[jy]= wy2*(wx2*zab[1,jy] + wx3*zab[2,jy])+ wy3*(wx2*zab[1,jy+1] + wx3*zab[2,jy+1]) } if( ixm1 == 0 ) next # Definition of coordinates for interface to farbrc xx[1]= x4 xx[2]= x3 xx[3]= x3 xx[4]= x4 yy[1]= y4 yy[2]= y4 yy[3]= y3 yy[4]= y3 zz[1]= z[ix,iy] zz[2]= z[ix-1,iy] zz[3]= z33 zz[4]= z43 zzx[1]= zx[2] zzy[1]= zy[2] zzxy[1]= zxy[2] zzx[2]= zx34 zzy[2]= zy34 zzxy[2]= zxy34 zzx[3]= zx33 zzy[3]= zy33 zzxy[3]= zxy33 zzx[4]= zx[1] zzy[4]= zy[1] zzxy[4]= zxy[1] nsides= 3 if (ix == 2) nsides= 4 ifa= ifa +1 call farbrc (xx,yy,zz,zzx,zzy,zzxy,cn,icol,nc,mode,nsides) } call frbfcl(icol) } end # End of farb2d #--------------------------------------------------------------------------- # farbrc - Fill area for a bicubic function on a rectangle # # Description # See below and for farb2d. # # History # 16Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- # f ill ar ea for a b icubic function on a r e c tangle # * ** * * * # # t r i p algorithm a.preusser farb-e-2d version 2.1 10/1988 # # author: a. preusser # fritz-haber-institut der mpg # faradayweg 4-6 # d-1000 berlin 33 # # # this subroutine computes a bicubic function from the # values x,y,z,zx,zy,zxy given at the four vertices of # a rectangle, and plots contours for the z-values cn[i], # i=1,nc, using the colours icol[i]. # area filling, set by parameter mode, is an optional feature. # # input parameters - See below # ============================ # # output # ====== # is performed by calls to the subroutine usrplt # to be supplied by the user (an example for usrplt # is included). # # parameters of usrplt # subroutine usrplt (x,y,n,ncol,mode) # x,y real arrays of length n for # the coordinates of the polygon # to be plotted. # n number of points of polygon # ncol index defining the colour for # the area or the line # mode 1, line drawing # 0, fill area # # if a rectangle receives only one colour, # the area is not filled at once. # instead, a 'fill area buffer' is opened # or updated, until a rectangle with a # different colour is encountered. # therefore, if the next call to farbrc # is not for a right-hand-neighbor, # or if it is the last call, subroutine # frbfcl # must be called by the user # call frbfcl(icol) , # in order to clear the fill area buffer, # and to fill the area of the rectangle. # # denomination of the vertices and sides of the # rectangle # y # side(3) # vertex(2) * -------------------------------0-------- * vertex(1) # ( . ) # ( . ) # ( . ) # side(4) ( . ride ) side(2) # ( . ) # ( . ) # ( . ) # vertex(3) * ----------------------------------0----- * vertex(4) # side(1) x # # the sides are parallel to the cartesian x-y-system. # # ----------------------------------------------------------------- # end of user documentation # ----------------------------------------------------------------- # # some nomenclature # # station zero on a side # ride move from one station to another inside rect. # transfer move from one station to the next on side # trip sequence of rides and transfers # round trip successful trip that ended at its start # horror trip trip that does not find an end # journey sequence of trips starting from the same # type of stations (same value of istatz) # and having the same orientation. # # there may be three journeys. # the first two are counter-clockwise and # start at stations with istatz=0 and =2, # respectively. the third journey is carried # out only in case of numerical difficulties, # when areas are unfilled after the first two. # it starts at stations with istatz=1 or =0 and # is clockwise. # # #--------------------------------------------------------------------------- procedure farbrc(x,y,z,zx,zy,zxy,cn,icol,nc,mode,nsides) real x[4], y[4], z[4] # Coordinates of the vertices # in centimeters (or inches, if cmscal=2.54). # real arrays of length (4). # x(i),y(i),z(i), i=1,4 define the # position of vertex (i). # The sides of the rectangle must be parallel # to the x- and y-axis, # and the vertices must be ordered # counter-clockwise as is indicated below # (vertex 1 in the upper right corner). real zx[4], zy[4], zxy[4] # Derivatives of z at the vertices. # real arrays of length (4). real cn[nc] # z- values in ascending order # for the contour levels. # real array of length (nc) int icol[nc+1] # Indices of the colours or patterns # for the areas between the contour lines. # the values of icol are passed to usrplt. # icol[i] is used for the area, where # z > cn[i-1] and # z <= cn[i], # for i=2,nc. # areas, where z <= cn[1] # are filled with colour icol[1], # and areas, where z > icol[nc] # are filled with colour icol[nc+1]. int nc # Number of contour levels, nc <= 100 int mode # What mode to produce the graph: # 0 -> fill area only # 1 -> lines only # 2 -> fill area and lines # Declarations real a, b, c real cc[100] # scaled contour levels. real cl # Actual contour level. real cmscal # Variable for switching between cm and inch. real cn1, cnn real d, derno1 real dernor # sign of normal derivative for a ride (+1 or -1). real dernos, di2, di3 real dx[4], dy[4] # Differences of x and y. real e, eps, epsn real f1, f1f2, f2, fa, fb real hmin # length of shortest side of rectangle. real p0[4], p1[4], p2[4], p3[4] # Coefficients for the polynomials # on the 4 sides . # variables on sides 1 and 2 are # counter-clockwise, on sides 3 and 4 # clockwise. real p11, p12, p13, p21, p22, p23, p31, p32, p33 # coeff. of polynomials used together # with p0[i]...p3[i], i=1 and i=4, for the # representation inside the rectangle. real pi # Value of PI, 3.141... real poserr # permitted position error. real q0[4], q1[4], q2[4] # Coefficients for the derivatives # of the polynomials on the 4 sides. real r0[4], r1[4] # Coeff. for the second derivatives. real sa[4], se[4] # Values of variables at start and end of sides. real sacmin # Minimal distance of two points to be stored. real sdchek real sder[50, 4] # Derivative in direction of sides # (side direction= counter clockwise). real si[4], co[4] # Cosinus of direction for sides. real sigs[4] # Sign for sides (+1 or -1). real sir, cor # Cosinus of direction normal to curve direction. real sl[4] # side lengths. real sl1[4] # 1./side length. real sl12[4] # sl1**2. real slmax # length of longest side. real t real t1[4, 4] # coordinates on sides at endpoints of intervals. real ta, tb real ts2[4] # working array for computing t1. real tzr[50, 4] # Coordinates for zeros (stations) on sides. real x0[50, 4], y0[50, 4] # x-y-coordinates of zeros on sides. real x3, y3 # coordinates x[3], y[3]. real xfabu[4], yfabu[4] # x,y coordinates of fill area buffer. real xpol[300], ypol[300] # x,y-coordinates for fill area polygon. real xstack[100], ystack[100] # x,y coordinates of stack # (coordinates of the first ride of a trip). real xx[2], yy[2] # local copies of x,y. real xx4f, yy4f # Coordinates for point p4f (preleminary position # of point p4). real z1[4, 4] # z-values at endpoints of intervals. real z3a3, zmaxt, zmint real zmin[4], zmax[4] # min. and max. values of z on sides. real zs real zsold # zs of last call. real zx3b3, zx4b3, zy3a3, zy4a3 real zz[4] # scaled z-values at vertices. int i, ic1, icn, icn1, ii int in[4] # number of intervals on sides. int iride, istart int istatz[50, 4] # status of zeros on sides # the initial status is = 0. # when a zero has served as start +1 is added. # when a zero has served as end +2 is added. int it int j, jj, jin, jjz, jn, journy, jpol, jpoll1, jr, js, jsa, jsa1st, jsa2 int jsar[6], jzar[6] # side and zero for start of rides. int jse, jse1st int jser[6], jzer[6] # side and zero for end of rides. int jsides int jstopr[6] # stop modes of rides. int jz, jza, jza1st, jza2, jze, jze1st, jzn, jzr int k, kcl, kcll int kk # Index of function to be evaluated by frbeva. int kride # counts the calls to frbrid. int kse # Actual side index. int kvert int maxpol # Maximum number of points for a trip. int maxrid # Maximum number of rides for a trip. int maxsta # Maximum number of points to be stored for a ride. int nc1, ncdif, ncdifa int nclzr[50, 4] # Contour level for zeros. int ncmax # Maximum number of contour levels. int ncmaxs # Maximum number of contours crossing a # rectangle side. int ncolbu # Color of fill area buffer. int ncol, ncl1 int ncpmax # Maximum number of points to be computed for a ride int ndir int ndir3 # 1, for counter-clockwise trip # -1, for clockwise trip int ndirs, ndirv int nfabu # Fill area buffer state. Possible values: # 0 -> fill area buffer closed. # 1 -> fill area buffer open. int nfar, ni, njour int np, np1, np1st, np2, npol1 int npp # Accumulated number of points for a sequence of # rectangles. int nprec # number of points for the rectangle. int nse, nside, nsides, nstop int nz[4] # Number of zeros on sides. # Function declarations. real frbeva(), frbzer() common /frbcob/ nfabu,ncolbu,xfabu,yfabu common /frbcoc/ sacmin,cmscal,maxpol,ncpmax,maxsta,ncmaxs,npp,pi, maxrid, sigs, ncmax # /frbcof/ contains variables which are passed to function # frbeva as parameters common /frbcof/ kk,kse,xx4f,yy4f,sir,cor,cl common /frbcop/ p0,p1,p2,p3, q0,q1,q2, r0,r1, p11,p12,p13,p21,p22,p23,p31,p32,p33 # frbcrd contains variables that are passed to frbrid # or that are retained for the next call to farbrc (nside=3) common /frbcrd/ x0,y0,nclzr,sder,tzr, nz,si,co,sa,se,dx,dy,sl, hmin,slmax,kride,nprec,poserr,dernor,ndir3, zmax,zmin,istatz,x3,y3,zsold #save it data xstack,ystack /200*0./ data it /0/, zs/0./ begin # initialisation for first rectangle it= it + 1 # set installation parameters # if (it == 1) { # cmscal is something for the Calcomp library to help in scaling. For # IRAF GIO though, which is capable of arbitrary coordinates, this should # always be left at 1. cmscal= 1. sacmin= 0.02/cmscal maxpol= 300 ncmaxs= 50 ncmax= 100 maxsta= 100 maxrid= 6 ncpmax= maxsta*10 npp= 0 nfabu= 0 pi= 4.*atan(1.) sigs[1]= 1. sigs[2]= 1. sigs[3]= -1. sigs[4]= -1. } # Check number of contour levels if (nc > ncmax) call error( FATAL_ERROR, "farbrc: number of contours to large" ) # nfar = number of fill area calls nfar= 0 # kride = number of calls to frbrid kride= 0 # nprec = number of curve points computed for rectangle nprec= 0 # check contour values for monotony, # if contour line passes through data point on # side 4, set nside to 4 (no copy from last rectangle) # scale contour levels nside= nsides if (cn[1] == z[2] || cn[1] == z[3]) nside= 4 zsold= zs zs= (z[1]+z[2]+z[3]+z[4])/4. cc[1]= cn[1] - zs do kcl=2,nc { cc[kcl]= cn[kcl] - zs if (cn[kcl] == z[2] || cn[kcl] == z[3]) nside= 4 if (cn[kcl]-cn[kcl-1] <= 0.) call error( FATAL_ERROR, "farbrc: contours are out of order" ) } # some basic geometry for the rectangle sl[1]= x[4]-x[3] sl[2]= y[1]-y[4] sl[3]= x[1]-x[2] sl[4]= y[2]-y[3] x3= x[3] y3= y[3] do j=1,4 { zz[j]= z[j] - zs np1= mod(j+1,4) + 1 np2= mod(j+2,4) + 1 dx[j]= x[np2] - x[np1] dy[j]= y[np2] - y[np1] sl1[j]= 1./sl[j] sl12[j]= sl1[j]*sl1[j] co[j]= dx[j]/sl[j] si[j]= dy[j]/sl[j] sa[j]= 0. if (j > 2) sa[j]= sl[j] se[j]= sl[j] if (j > 2) se[j]= 0. } slmax= amax1(sl[1],sl[2]) di2= - (dy[1]*co[2] - dx[1]*si[2]) di3= abs(dy[2]*co[3] - dx[2]*si[3]) # check coordinates of vertices if (y[1] != y[2] || x[2] >= x[1]) call error( FATAL_ERROR, "farbrc: rectangle vertices not ordered counter-clockwise") # check if vertices are numbered counter-clockwise if (di2 < 0.) call error( FATAL_ERROR, "farbrc: rectangle vertices not ordered counter-clocwise" ) # = shortest side length hmin= amin1(di2,di3] # check hmin if (hmin < 0.01/cmscal || hmin > 100./cmscal) call error( FATAL_ERROR, "farbrc: differences in x,y coordinates too small or large to scale" ) if (hmin == 0.) return # = permitted position error poserr= amin1(1.e-03/cmscal,1.e-03*hmin*hmin/slmax) # copy information for side 4 if (nside != 4) { jzn= nz[2] nz[4]= jzn zmin[4]= zmin[2] + zsold - zs zmax[4]= zmax[2] + zsold - zs if (jzn != 0) { jjz= jzn do jzr=1,jzn { x0[jjz,4]= 0. y0[jjz,4]= y0[jzr,2] nclzr[jjz,4]= nclzr[jzr,2] sder[jjz,4]= -sder[jzr,2] tzr[jjz,4]= tzr[jzr,2] istatz[jjz,4]= 0 jjz= jjz - 1 } } } # Compute coefficients for polynomials along sides z3a3= (zz[4]- zz[3])*sl1[1] p0[1]= zz[3] p1[1]= zx[3] p2[1]= (2.0*(z3a3-zx[3])+z3a3-zx[4])*sl1[1] p3[1]= (-2.*z3a3+zx[4]+zx[3])*sl12[1] z3a3= (zz[2]- zz[3])*sl1[4] p0[4]= zz[3] p1[4]= zy[3] p2[4]= (2.*(z3a3-zy[3])+z3a3-zy[2])*sl1[4] p3[4]= (-2.*z3a3+zy[2]+zy[3])*sl12[4] z3a3= (zz[1]- zz[2])*sl1[3] p0[3]= zz[2] p1[3]= zx[2] p2[3]= (2.*(z3a3-zx[2]) + z3a3 - zx[1])*sl1[3] p3(3]= (-2.*z3a3+zx[1]+zx[2])*sl12[3] z3a3= (zz[1]- zz[4])*sl1[2] p0[2]= zz[4] p1[2]= zy[4] p2[2]= (2.*(z3a3-zy[4]) + z3a3 - zy[1])*sl1[2] p3[2]= (-2.*z3a3+zy[1]+zy[4])*sl12[2] # determine polynomial coeff for derivatives along sides do j= 1,4 { q0[j]= p1[j] q1[j]= 2.*p2[j] q2[j]= 3.*p3[j] r0[j]= q1[j] r1[j]= 2.*q2[j] } # set contour level to be passed to frbeva cl= 0. # find points on sides of rectangle, # where first derivative is zero # # loop over sides do jsa=1,nside { kse= jsa # set initial endpoints of intervals t1[1,kse]= sa[kse] t1[2,kse]= se[kse] i= 2 # loop over derivatives do k=3,4 { # set function to be evaluated by frbeva kk= k ts2[1]= t1[1,kse] ii= 2 tb= t1[1,kse] f2= frbeva(tb) # loop over endpoints of intervals do j=2,i { ta= tb f1= f2 tb= t1[j,kse] f2= frbeva(tb) if (f1*f2 > 0.) next if (f1 == 0. && f2 == 0.) next ts2[ii]= frbzer(ta,tb,f1,f2,poserr) ii= ii + 1 } ts2[ii]= t1[i,kse] do j=1,ii t1[j,kse]= ts2[j] i= ii } # in[kse]= number of intervals # (e.g. if in[kse]=1, there is no point for which 1st der.=0) i= i-1 in[kse]= i # compute maxima and minima for each side np1= mod(kse+1,4) + 1 np2= mod(kse+2,4) + 1 zmax[kse]= amax1(zz[np1],zz[np2]) zmin[kse]= amin1(zz[np1],zz[np2]) z1[1,kse]= zz[np1] z1[i+1,kse]= zz[np2] if (i != 1) { kk= 1 do j=2,i { z1[j,kse]= frbeva(t1[j,kse]) if (z1[j,kse] > zmax[kse]) zmax[kse]= z1[j,kse] if (z1[j,kse] < zmin[kse]) zmin[kse]= z1[j,kse] } } } # check, if rectangle has one colour, # because the minimum is over the max. contour level, # or maximum under min. contour level zmaxt= amax1(zmax[1],zmax[2],zmax[3],zmax[4]) zmint= amin1(zmin[1],zmin[2],zmin[3],zmin[4]) cn1= cc[1] cnn= cc[nc] if (cnn < zmint || cn1 >= zmaxt ) { ncol =1 if (cnn < zmint) ncol= nc+1 label_25_ # rectangle has one colour only if (mode != 1) { if (nfabu == 1 && ncol == ncolbu) call frbfup (x,y) else call frbfop(x,y,icol,ncol) } nz[2]= 0 return } # find min. and max. contour level for rectangle do kcl=1,nc { ic1= kcl if (cc[kcl] >= zmint) break } icn= nc + 1 do kcl=1,nc { icn= icn - 1 if (cc[icn] <= zmaxt) break } # ic1= first contour level # icn= last contour level icn1= icn - ic1 + 1 ncol= icn +1 if (cc[icn] == zmaxt) ncol= icn if (icn1 == 0) goto label_25_ # compute zeros on sides for all contour levels # in counter-clockwise order do jsa=1,nside { kse= jsa jn= 0 if (cc[icn] >= zmin[kse] && cc[ic1] <= zmax[kse] ) { ni= in[kse] do jin=1,ni { ndir= sign(1.,z1[jin+1,kse]-z1[jin,kse]) kcl= ic1 - 1 if (ndir == -1) kcl= icn + 1 do kcll=1,icn1 { kcl= kcl + ndir cl= cc[kcl] f1= z1[jin,kse] - cl f2= z1[jin+1,kse] - cl f1f2= f1*f2 if (f1f2 > 0.) next # special situations if (f1f2 == 0.) { # if () then contourline = side kse if (ni == 1 && f1 == 0. && f2 == 0.) goto label_67_ # if () then line passes through a vertex at start of side # this case is handled on previous side if (f1 == 0. && abs(t1[jin,kse]-t1[1,kse]) <= poserr*3.) next # if () then line passes through vertex at end of side if (f2 == 0. && abs(t1[jin+1,kse]-t1[ni+1,kse]) <= poserr*3.) goto label_68_ goto label_69_ # contour line = side jsa label_67_ kvert= mod(jsa+1,4) + 1 xx[1]= x[kvert] yy[1]= y[kvert] kvert= mod(jsa+2,4) + 1 xx[2]= x[kvert] yy[2]= y[kvert] if (mode > 0) call usrplt(xx,yy,2,icol[kcl],1) # call usrplt(xx,yy,2,kcl,1) goto label_85_ label_68_ # line passes through data point # (at end of side) # inhibit multiple zero at vertex if ( jn != 0 && kcl == nclzr[jn,kse] ) if (abs(t1[ni+1,kse]-tzr[jn,kse]) <= poserr*3.) next # compute value on side jsa kk= 1 eps= 0.01*sl[kse] ta= t1(ni+1,kse]-eps*sigs[kse] fa= frbeva(ta) # compute value on next side nse= mod(kse,4) + 1 kse= nse epsn= 0.01*sl[kse] tb= t1[1,nse]+ epsn*sigs[nse] fb= frbeva(tb) kse= jsa # if () then contour line is degenerated to a point if (fa*fb > 0.) next # contour line starts from vertex into rectangle jn= jn+ 1 if (jn > ncmaxs) call error( FATAL_ERROR, "farbrc: too many contours crossing a rectangle boundary" ) tzr[jn,kse]= t1[ni+1,kse] nclzr[jn,kse]= kcl sder[jn,kse]= -fa/eps next } # compute zero on side (station) label_69_ jn= jn+1 if (jn > ncmaxs) call error( FATAL_ERROR, "farbrc: too many contours crossing a rectangle boundary" ) kk= 1 tzr[jn,kse]= frbzer(t1[jin,kse],t1[jin+1,kse],f1,f2, poserr) # compute derivative at zero kk= 4 sder[jn,kse]= frbeva(tzr[jn,kse])*sigs[kse] # store index of contour level nclzr[jn,kse]= kcl # check sign of der., if same level if (jn < 2) next if (kcl != nclzr[jn-1,kse]) next if (abs(tzr[jn,kse]-tzr[jn-1,kse]) > poserr*3.) next if (sder[jn,kse]*sder[jn-1,kse] < 0.) next # if sign is wrong or =0, compute der. by differences kk= 1 eps= 0.01*sl[kse] ta= tzr[jn-1,kse] - eps*sigs[kse] sder[jn-1,kse]= -frbeva(ta)/eps tb= tzr[jn,kse] + eps*sigs[kse] sder[jn,kse]= frbeva(tb]/eps } } } label_85_ # nz = number of zeros on side kse nz[kse]= jn } # every ride should have start and end if (nz[1]+nz[2]+nz[3]+nz[4] < 2) goto label_25_ # clear fill area buffer call frbfcl (icol) # compute x0,y0 for each zero (relative to x[3],y[3]), # set status of all zeros to 0 do jsa=1,nside { jn= nz[jsa] if (jn == 0) next np1= mod(jsa+1,4) + 1 do jza=1,jn { istatz[jza,jsa]= 0 t= tzr[jza,jsa] x0[jza,jsa]= x[np1] -x3 + sa[jsa]*co[jsa] +t*abs(co[jsa]) y0[jza,jsa]= y[np1] -y3 + sa[jsa]*si[jsa] +t*abs(si[jsa]) } } # compute coefficients for representation inside rectangle zx3b3= (zx[2]-zx[3])*sl1[4] zx4b3= (zx[1]-zx[4])*sl1[4] zy3a3= (zy[4]-zy[3])*sl1[1] zy4a3= (zy[1]-zy[2])*sl1[1] a= (zz[1]-zz[4]-zz[2]+zz[3])*sl1[4]*sl1[1] - zx3b3 - zy3a3 + zxy[3] b= zx4b3 - zx3b3 - zxy[4] + zxy[3] c= zy4a3 - zy3a3 - zxy[2] + zxy[3] d= zxy[1] - zxy[4] - zxy[2] + zxy[3] e= a+a-b-c p11= zxy[3] p12= (2.*(zx3b3-zxy[3])+zx3b3-zxy[2])*sl1[4] p13= (-2.*zx3b3+zxy[2]+zxy[3])*sl12[4] p21= (2.*(zy3a3-zxy[3])+zy3a3-zxy[4])*sl1[1] p22= (3.*(a+e)+d)*sl1[1]*sl1[4] p23= (-3.*e-b-d)*sl1[1]*sl12[4] p31= (-2.*zy3a3+zxy[4]+zxy[3])*sl12[1] p32= (-3.*e-c-d)*sl1[4]*sl12[1] p33= (d+e+e)*sl12[1]*sl12[4] # initialize stack (for first ride of a trip) jse1st= 0 jze1st= 0 jsa1st= 0 jza1st= 0 np1st= 0 dernos= 0 # start 'trips' using the zeros on the sides as 'stations' # # set parameters for first journey # start at unused stations (istatz =0) istart= 0 # ndir3= 1 means counter-clockwise trip ndir3= 1 ndirv= 2 ndirs= 0 # loop over journeys njour= 3 if (mode == 1) njour=1 do journy= 1,njour { # set parameters for journy=3 if (journy == 3) { istart= 1 ndir3= -1 ndirv= 1 ndirs= 2 } # loop over sides, jsa= side index do jsa=1,4 { jzn= nz[jsa] if (jzn != 0) { # loop over zeros, jza= starting zero do jza= 1,jzn { # for third journey, start also at istatz=0 if (journy != 3 || istatz[jza,jsa] != 0) # if () then this station will not serve as start if (istatz[jza,jsa] != istart) next # start trip from station side jsa, zero jza # if journy != 1 check stack first if (journy != 1 && jsa == jse1st && jza == jze1st) { # coordinates for next ride are in stack jj= np1st do j=1,np1st { xpol[j]= xstack[jj] ypol[j]= ystack[jj] jj= jj-1 } npol1= np1st jsa2 = jsa1st jza2 = jza1st nstop= 0 # set normal deriv. and color dernor= -dernos derno1= dernor ncl1= nclzr[jza,jsa] ncol= amin1(-dernor*ndir3+1.,1.) + ncl1 } else { # set normal derivative and colour for first ride if (sder[jza,jsa] == 0.) next dernor= sign(1.,sder[jza,jsa]) derno1= dernor ncl1= nclzr[jza,jsa] ncol= amin1(-dernor*ndir3+1.,1.) + ncl1 # ncol= colour for trip # the first line of the last statement is 1 or 0, # depending on the sign of sder # (derivative in direction of side) # and ndir3 # ndir3= 1, for journy=1,2 (counter-clockwise trip) # ndir3= -1, for journy=3 (clockwise trip) # call frbrid(jsa,jza,cc,nc,xpol,ypol,maxpol,jsa2,jza2,npol1, nstop) # first ride ended on station side jsa2, zero jza2 # # call for line drawing if (mode == 1 ) call usrplt(xpol,ypol,npol1,icol[ncl1],1) # call usrplt(xpol,ypol,npol1,ncl1,1) } # book keeping for start and end of ride jpol= npol1 jpoll1= 1 np= npol1 iride= 1 jser[1]= jsa2 jzer[1]= jza2 jsar[1]= jsa jzar[1]= jza jstopr[1]= nstop jsides= 0 if (nstop == 2) next # find next station for continuation of trip (=transfer) label_36_ jza2 = jza2 + ndir3 if (jza2 > nz[jsa2] || jza2 <= 0) { # take next side, add vertex to polygon # set jza2 to first or last zero of the new side repeat { jsides= jsides + 1 if (jsides == 5) next jpol= jpol + 1 kvert= mod(jsa2+ndirv,4) +1 xpol[jpol]= x[kvert] ypol[jpol]= y[kvert] jsa2= mod(jsa2+ndirs,4) + 1 } until (nz[jsa2] != 0) jza2= 1 if (ndir3 == -1) jza2= nz[jsa2] } # check for regular end of trip if (jsa2 != jsa || jza2 != jza) { # check difference of contour levels # between jsa,jza and jsa2,jza2, # check sign of derivative of station jza2,jsa2, # set new dernor # ncdif= nclzr[jza2,jsa2] - ncl1 ncdifa= abs(ncdif) nc1= -ncdifa if (nc1 == 0) nc1= 1 dernor= derno1*nc1 sdchek= sder[jza2,jsa2]*dernor if (sdchek < 0. || ncdifa > 1) { # do not stop at vertex if (tzr[jza2,jsa2] < poserr*3. || tzr[jza2,jsa2] > sl[jsa2]-poserr*3.) goto label_36_ # stop trip in all other cases next } # start new ride from side jsa2, zero jza2 # # check stack first if (jsa2 == jse1st && jza2 == jze1st) { # coordinates for next ride are in stack # the following replaces a call to frbrid if (jpol+np1st > maxpol) call error( FATAL_ERROR, "farbrc: Overflow of working storage xpol, ypol" ) jj= np1st do j=1,np1st { xpol[jpol+j]= xstack[jj] ypol[jpol+j]= ystack[jj] jj= jj-1 } jse= jsa1st jze= jza1st np= np1st nstop= 0 } else { call frbrid (jsa2,jza2,cc,nc, xpol[jpol+1],ypol[jpol+1],maxpol-jpol,jse,jze,np,nstop) # ride ended at station side jse, zero jze # # call for line drawing if (mode == 1 && istatz[jza2,jsa2] == 0 && np >= 2) call usrplt(xpol[jpol+1],ypol[jpol+1], np,icol[nclzr[jza2,jsa2]],1) # call usrplt(xpol[jpol+1],ypol[jpol+1], # np,nclzr[jza2,jsa2],1) } # book keeping for continuation ride if (nstop == 2) np= 0 jpoll1= jpol + 1 jpol= jpol + np iride= iride + 1 jsar[iride]= jsa2 jzar[iride]= jza2 jser[iride]= jse jzer[iride]= jze jstopr[iride]= nstop jsa2= jse jza2= jze if (iride+1 > maxrid) next # continue trip goto label_36_ } # polygon for fill area is complete # (successfull round trip) # # write coordinates of first ride to stack if ( npol1 <= maxsta && jstopr[1] != 1 ) { do j=1,npol1 { xstack[j]= xpol[j] ystack[j]= ypol[j] } jsa1st= jsa jza1st= jza jse1st= jser[1] jze1st= jzer[1] dernos= derno1 np1st= npol1 } # call for fill area if (jpol < 3) next nfar= nfar + 1 if (mod(mode,2) == 0) call usrplt(xpol,ypol,jpol,icol[ncol],0) # set flags for start and end of rides do jr= 1,iride { js= jsar[jr] jz= jzar[jr] if (mod(istatz[jz,js],2) == 0) istatz[jz,js]= istatz[jz,js] + 1 if (jstopr[jr] > 0) next js= jser[jr] jz= jzer[jr] if (istatz[jz,js] < 2) istatz[jz,js]= istatz[jz,js] + 2 } # draw line if mode=2 js= jsar[iride] jz= jzar[iride] if (mode == 2 && np > 1) call usrplt(xpol[jpoll1],ypol[jpoll1],np,icol[nclzr[jz,js]],1) # call usrplt(xpol[jpoll1],ypol[jpoll1],np,nclzr[jz,js],1) } # end of loop over zeros } } # end of loop over sides # # set istart for journy=2 istart= 2 } npp= npp + nprec # treat case when whole rectangle has to be filled # because there was no successfull trip if (nfar == 0 && mode != 1) { ncol= icn + 1 if (cc[icn] == zmaxt) ncol= icn goto label_25_ } end #--------------------------------------------------------------------------- # End of farbrc #--------------------------------------------------------------------------- # frbfcl - Clear fill area buffer # # Description # See above in farb2d. # # History # 17Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- procedure frbfcl(icol) int icol[*] # Declarations real xfabu[4], yfabu[4] int ncolbu, nfabu common /frbcob/ nfabu,ncolbu,xfabu,yfabu begin # If the fill area isn't closed, close it and write it out, else do nothing. if (nfabu == 1) { # declare fill area buffer closed nfabu= 0 # call fill area routine call usrplt(xfabu,yfabu,4,icol[ncolbu],0) } end #--------------------------------------------------------------------------- # End of frbfcl #--------------------------------------------------------------------------- # frbrid - trace contour from side jsa to side jsa2 (ride from jsa,jza to jsa2,jza2) # Description # See above in farb2d and below. # # Author # t r i p algorithm a.preusser farb-e-2d version 2.1 10/1988 # # # author: a. preusser # fritz-haber-institut der mpg # faradayweg 4-6 # d-1000 berlin 33 # # History # 17Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- procedure frbrid (jsa,jza,cn,nc,xpol,ypol,maxpol,jsa2,jza2,np, nstop) int jsa # I: side index contour starts from int jza # I: zero index contour starts from real cn[nc] # I: contour levels real xpol[maxpol], ypol[maxpol] # O: x-y-coordinates of the points of a ride int maxpol # I: maximum number of points in xpol,ypol int jsa2 # O: side index where contour ends int jza2 # O: zero (station) index where contour ends int np # O: number of points stored to xpol,ypol int nstop # O: Where the ride ended: values are: # 0 -> ride ended at station # 1 -> ride ended on side, no station found # 2 -> ride ended inside rectangle # If nstop == 1, jsa2,jza2 indicate the # previous station on the ending side. # (note positive or negative sense of trip # indicated by ndir3). jza2 may be zero. # # If nstop == 2, jsa2= jsa; jza2= jza . # Declarations real cl, cmscal, co[4], cor real dernor, ds, ds12, ds13, ds23, ds30, ds34, dsmax, dsmin, dx[4], dx12, dx30, dxds, dxds1, dy[4], dy12, dy30, dyds, dyds1 real eps real f1, f2, fstep real hmin real ogr real pi, pl0, pl1, pl2, poserr real r, r30, r30min, racc, racmin, rma, rmax real s30max, sa[4], sacc, sacmin, sder[50,4], se[4], si[4], sigs[4], sir, sl[4], slmax, soll1, sq real thetac, thetas[4], thetsc, ttzr, tzr[50,4] real ugr real x0[50,4], x3, xx1, xx2, xx3, xx4, xx4f real y0[50,4], y3, yy1, yy2, yy3, yy4, yy4f real zmax[4], zmin[4], zsold int istatz[50,4] int j, jese, jj, jn, jp, jpp, jse int kcl, kk, kride, kse int maxrid, maxsta int n, nc, nclzr[50,4], ncmax, ncmaxs, ncp, ncpmax, ndir3, nost, npmax, npp, nprec, nsad, nse, nz[4] # Function declarations. real frbeva(), cos(), sin(), frbzer() common /frbcoc/ sacmin,cmscal,npmax,ncpmax,maxsta,ncmaxs,npp,pi, maxrid, sigs, ncmax common /frbcof/ kk,kse,xx4f,yy4f,sir,cor,cl common /frbcrd/ x0,y0,nclzr,sder,tzr, nz,si,co,sa,se,dx,dy,sl, hmin,slmax,kride,nprec,poserr,dernor,ndir3, zmax,zmin,istatz,x3,y3,zsold common /frbcor/ rma,rmax,dsmax,dsmin,fstep, thetas,racmin #save /frbcor/ begin kse= jsa kride= kride +1 nsad= 0 nstop= 0 # Initialization, if first ride of a rectangle if (kride == 1) { rmax = amin1(0.010/cmscal,hmin*0.010) # = distance normal to curve direction within which a zero must be found dsmax = hmin*0.2 # = maximum step size dsmin = amin1(rmax*0.03,poserr*8.) # = minimum step size fstep = amin1(rmax*8.,dsmax) # = starting step size racmin= sacmin*0.1 # a point is stored only if the accumulated r's (racc) # (change in direction) have reached racmin # set direction of sides thetas[1]= 0. thetas[2]= pi*0.5 thetas[3]= pi thetas[4]= pi*1.5 } # define contour level kcl= nclzr[jza,kse] cl= cn[kcl] # estimate starting direction # # compute f1 on side nse= mod(kse,4) + 1 eps= 0.01*sl[nse] f1= sder[jza,kse]*eps # compute f2 normal to side kk= 2 xx4f= x0[jza,kse] yy4f= y0[jza,kse] sir= co[kse] cor= - si[kse] f2= frbeva(eps) # compute angle for starting direction if (f2 == 0.) { # if (f1 == 0.) goto label_69_ # NOTE: The above was commented out in the orignal source also: jde thetsc= pi*0.5 } else { thetsc= atan (-f1/f2) if (thetsc < 0.) thetsc= thetsc + pi if (thetsc == 0.) thetsc= (ndir3-1)*pi*0.5 } thetac= thetas[kse] + thetsc # compute points # # store first point, jp= number of points stored for this contour line jp= 1 xpol[jp]= x0[jza,kse] ypol[jp]= y0[jza,kse] # initialize tracing label_49_ r= 0. ds= fstep dx12= ds*cos(thetac) dy12= ds*sin(thetac) xx3= xx4f yy3= yy4f xx2= xx3 - dx12 yy2= yy3 - dy12 xx1= xx2 - dx12 yy1= yy2 - dy12 ds23= ds ds12= ds # ds01= ds # points p1,p2,p3,p4 with coordinates xx1...xx4, yy1...yy4 # are referred to as *queue*. ds12...ds34 are the distances # between points in the queue. # xx4f,yy4f are preleminary coordinates for the next point p4 # which will be computed by the regula falsi (frbzer). # for a derivation of the formulas for pl0...pl2 see # preusser,a. computing area filling contours for surfaces # defined by piecewise polynomials. # computer aided geometric design 3, # (1986), p. 267-279 # there is also an explanation for the following part of # the algorithm. # = accumulated distances to last point stored sacc= 0. # = accumulated r # a point is only stored to xpol,ypol if sacc >= sacmin # and racc >= racmin racc= racmin # = number of points computed ncp= 0 # = number of steps normal to curve nost= 0 rma= rmax # compute new point for contour line # # compute curve direction label_15_ ds13= ds23 + ds12 pl0= ds23/(ds12*ds13) pl1= -ds13/(ds12*ds23) pl2= (ds13+ds23)/(ds13*ds23) dxds= pl0*xx1 + pl1*xx2 + pl2*xx3 dyds= pl0*yy1 + pl1*yy2 + pl2*yy3 sq= sqrt(dxds*dxds+dyds*dyds) dxds= dxds/sq dyds= dyds/sq cor= -dyds sir= dxds # # search for two points with opposite sign rma= sign (rmax,rma) repeat { xx4f= xx3 + dxds*ds yy4f= yy3 + dyds*ds f1= frbeva(0.) rma= sign(rma,dernor*f1) f2= frbeva(rma) if (f1*f2 <= 0.) goto label_16_ label_56_ if (ds*0.5 < dsmin) break # divide stepsize in curve direction by 2. ds= ds*0.5 } # divide stepsize normal to curve by 2. repeat { nost= nost + 1 rma= rma*0.5 if (abs(rma) <= poserr) break f2= frbeva(rma) if (f1*f2 <= 0.) goto label_16_ } # saddle point # # set new direction nsad= nsad + 1 if (nsad > 1) goto label_69_ dxds1= dxds dyds1= dyds thetac= atan2(dyds,dxds) + pi*0.5 # store saddle point xx4f= xx3 yy4f= yy3 jpp= jp +1 if (jpp > maxpol) goto label_40_ jp= jpp xpol(jp]= xx3 ypol[jp]= yy3 goto label_49_ # find zero for new point label_16_ r= frbzer(0.,rma,f1,f2,poserr) ncp= ncp + 1 if (ncp > ncpmax) goto label_69_ ds34= sqrt(ds*ds + r*r) xx4= xx4f + cor*r yy4= yy4f + sir*r # check if point is outside the rectangle if (yy4 < 0. ) { jsa2=1 goto label_17_ } else if (xx4 > dx[1] ) { jsa2=2 goto label_17_ } else if (yy4 > dy[2] ) { jsa2= 3 goto label_17_ } else if (xx4 < 0. ) { jsa2= 4 goto label_17_ } else { # point is inside # # store point to xpol,ypol sacc= sacc + ds34 racc= racc + abs(r) if (sacc >= sacmin && racc >= racmin) { jpp= jp + 1 if (jpp > maxpol) goto label_40_ jp= jpp xpol[jp]= xx4 ypol[jp]= yy4 sacc= 0. racc= 0. } # update queue # ds01= ds12 ds12= ds23 ds23= ds34 xx1= xx2 yy1= yy2 xx2= xx3 yy2= yy3 xx3= xx4 yy3= yy4 # set new step size soll1= 2. if (abs(r) > poserr) soll1= abs(rma*0.8/r) if (soll1 > 1) ds= amin1(dsmax,ds*sqrt(soll1)) # goto label_15_ } # tracing stopped label_69_ nstop= 1 jsa2= 1 dxds= dxds1 dyds= dyds1 # point is outside # # search for corresponding zero on sides # start with side jsa2 label_17_ r30min= 99999. s30max= dsmax jse= jsa2 do n=1,2 { do jese=1,4 { jj= nz[jse] if (jj != 0) do j=1,jj { # level check if (nclzr[j,jse] != nclzr[jza,kse]) next # derivative check if (sder[j,jse]*dernor > 0.) next # r-check dx30= x0[j,jse] - xx3 dy30= y0[j,jse] - yy3 r30= (dy30*dxds - dx30*dyds) if (r30*r < 0 && n == 1) next r30= abs(r30) if (r30 >= r30min) next # s-check ds30= dx30*dxds + dy30*dyds if (ds30 > s30max || ds30 < -dsmin*2.) next jza2= j jsa2= jse r30min= r30 } jse= mod(jse,4) + 1 } # for n=2, do not check sign of r30 if (r30min < rmax) { nstop = 0 goto label_85_ } } # no acceptable zero on all three sides # # reduce step size if (nstop != 1) goto label_56_ label_85_ if (jp+1 > maxpol) goto label_40_ if (nstop != 1) { jp= jp+1 xpol[jp]= x0[jza2,jsa2] ypol[jp]= y0[jza2,jsa2] } else { # error handling if (jp <= 2) goto label_39_ # check if last point is near boundary if (abs(yy3) < sacmin) { jsa2= 1 ttzr= xx3 } else if (abs(xx3-dx[1]) < sacmin) { jsa2= 2 ttzr= yy3 } else if (abs(yy3-dy[2]) < sacmin) { jsa2= 3 ttzr= xx3 } else if (abs(xx3) < sacmin) { jsa2= 4 ttzr= yy3 } else goto label_39_ # tracing stopped near boundary # nstop= 1, jsa2 side number, # jza2 zero index of last or next zero, or =0 # depending on ndir3. # find jza2 ugr= sa[jsa2] jza2= 0 jn= nz[jsa2] if (jn != 0) { do jj=1,jn { jza2= jza2 + 1 ogr= tzr[jza2,jsa2] if (ttzr >= ugr && ttzr <= ogr || ttzr >= ogr && ttzr <= ugr) { jza2 = jza2 - 1 break } ugr= ogr } jza2= jza2 + 1 if (ndir3 == 1) jza2= jza2-1 } } # add coordinates of lower left vertex # normal return label_20_ np= jp do jp=1,np { xpol[jp]= xpol[jp] + x3 ypol[jp]= ypol[jp] + y3 } nprec= nprec + np return # tracing stopped within rectangle label_39_ nstop= 2 jsa2= kse jza2= jza goto label_20_ label_40_ call eprintf( "frbrid: warning: overflow of working storage xpol, ypol" ) goto label_39_ end #--------------------------------------------------------------------------- # End of frbrid #--------------------------------------------------------------------------- # frbfop - open fill area buffer # # Description # See farb2d. # # History # 17Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- procedure frbfop(x,y,icol,ncol) real x[4], y[4] int ncol, icol[ncol] # Declarations real xfabu[4], yfabu[4] int j, ncolbu, nfabu common /frbcob/ nfabu,ncolbu,xfabu,yfabu begin # if already open, clear buffer first if (nfabu == 1) call frbfcl (icol) # fill buffer do j=1,4 { xfabu(j)= x(j) yfabu(j)= y(j) } # set colour of fill area buffer ncolbu= ncol # declare fill area buffer open nfabu= 1 end #--------------------------------------------------------------------------- # End of frbfop #--------------------------------------------------------------------------- # frbfup - update fill area buffer # # Description # See farb2d. # # History # 17Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- procedure frbfup(x,y) real x[4], y[4] # Declarations real xfabu[4], yfabu[4] int ncolbu, nfabu common /frbcob/ nfabu,ncolbu,xfabu,yfabu begin xfabu[4]= x[4] yfabu[4]= y[4] xfabu[1]= x[1] yfabu[1]= y[1] end #--------------------------------------------------------------------------- # End of frbfup #--------------------------------------------------------------------------- # frbzer - Computer zero between limits. # # Description # The method is a combination of the regula falsi # and the midpoint method # # It is a modified version of the vim- (control data # user group) routine with catalog identification # c2bkyzero # written by loren p. meissner, 1965 # # Author # t r i p algorithm a.preusser farb-e-2d version 2.1 10/1988 # # History # 17Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- real procedure frbzer (ta,tb,f1,f2,er) real ta, tb # Limits real f1 # function value at ta real f2 # function value at tb # f1 and f2 must have opposite sign. # This must be checked before entry real er # Permitted error for solution frbzer # Declarations real a, b, c, e, fa, fb, fc, fg, fs, fy, g, h real s, y # Function declarations. real frbeva() begin a=ta b=tb fa=f1 fb=f2 c=a fc=fa s=c fs=fc repeat { h=0.5*(b+c) if(abs(h-b) <= er) break if (abs(fb) <= abs(fc)) { y=s fy=fs g=c fg=fc s=b fs=fb } else { y=b fy=fb g=b fg=fb s=c fs=fc } if (fy == fs) b=h else { e=(s*fy-y*fs)/(fy-fs) if (abs(e-s) <= er) e=s+sign(er,g-s) if ((e-h)*(s-e) < 0.0) b=h else b=e } fb=frbeva(b) if (fg*fb >= 0.0) { c=s fc=fs } else { c=g fc=fg } } return( h ) end #--------------------------------------------------------------------------- # End of frbzer #--------------------------------------------------------------------------- # frbeva - function evaluation # # Description # See farb2d. # # Author # t r i p algorithm a.preusser farb-e-2d version 2.1 10/1988 # # author : a. preusser # fritz-haber-institut # der max-planck-gesellschaft # faradayweg 4-6 # d-1000 berlin 33 # # History # 17Jan91 - Rewritten into SPP. Jonathan D. Eisenhamer, STScI. #--------------------------------------------------------------------------- real procedure frbeva(t) real t # Declarations. real answer real cl, cor real p0[4], p1[4], p2[4], p3[4] real p11, p12, p13, p21, p22, p23, p31, p32, p33 real q0[4], q1[4], q2[4] real r0[4], r1[4] real s0, s1, s2, s3, sir real xx4, xx4f real yy4, yy4f int kk # number of function to be evaluated. Possible values: # kk=1 original polynomial along side kse # kk=2 bivariate polynomial inside rectangle # kk=3 2nd derivative along side kse # kk=4 1st derivative along side kse int kse int kse # variabels in /frbcof/ are used as arguments # for an explanation see # subroutine farbrc common /frbcof/ kk,kse,xx4f,yy4f,sir,cor,cl common /frbcop/ p0,p1,p2,p3, q0,q1,q2, r0,r1, p11,p12,p13,p21,p22,p23,p31,p32,p33 begin switch( kk ) { case 2: xx4= xx4f + cor*t yy4= yy4f + sir*t s0= p0[1] + yy4*(p1[4]+yy4*(p2[4]+yy4*p3[4])) s1= p1[1] + yy4*(p11+yy4*(p12+yy4*p13)) s2= p2[1] + yy4*(p21+yy4*(p22+yy4*p23)) s3= p3[1] + yy4*(p31+yy4*(p32+yy4*p33)) answer= s0 + xx4*(s1+xx4*(s2+xx4*s3)) - cl case 3: answer= r0[kse] + t*r1[kse] case 4: answer= q0[kse] + t*(q1[kse]+t*q2[kse]) default: answer= p0[kse]+t*(p1[kse]+t*(p2[kse]+t*p3[kse])) - cl } return( answer ) end #--------------------------------------------------------------------------- # End of frbeva #---------------------------------------------------------------------------