program cisq(cisqp, xyin, output); (* cisq: circle to square Tom Schneider National Cancer Institute Laboratory of Mathematical Biology Frederick, Maryland 21701-1013 toms@ncifcrf.gov 1989 *) label 1; (* end of program *) const (* begin module version *) version = 1.43; (* of cisq.p 1989 December 19 origin 1989 November 24 *) (* end module version *) (* begin module describe.cisq *) (* name cisq: circle to square synopsis cisq(cisqp: in, xyin: out, output: out) files cisqp: parameters to control the program First line: lowest value of m, mlo. Second line: highest value of m, mhi. Third line: increment in the value of m, mstep. Fourth line: desired radius of a circle if m = 2, reffective. Fifth line: number of steps to take to move around 360 degrees. Sixth line: A factor by which to increase the value of theta, spinfactor. 1 gives a square, 1.5 gives a hexagon. xyin: input to the xyplo program. Curves that are close to integer values of n have the symbol m, others have the symbol r. This allow them to be distinguished by the graphics routines. output: messages to the user description Plot the equation |x|^m + |y|^m = |reffective|^m where reffective is the "effective" radius of the curve, |x| is the absolute value of x, and ^ means to raise to the mth power. This gives a line if m = 1, a circle if m = 2 and approaches a square as m -> infinity! The method for producing the curves is to re-express the equation in polar coordinates. One must be a bit careful to distinguish between the effective radius (reffective) and the current polar coordinate (r). After making this distinction we can write: x = r cos theta y = r sin theta and rearrange to solve for r, while keeping reffective fixed as it should be. Dividing the basic formula by r (>0) and converting to polar coordinates gives: (reffective/r)^m := / ((|cos(theta)|)^m + (|sin(theta)|)^m); To do this in Pascal, we have to use the form, a^m = exp(m*ln(a)). This gives: exp(m * ln(r/reffective)) := 1 / ( exp(m * ln(abs(cos(theta)))) + exp(m * ln(abs(sin(theta)))) ) where we have also introduced the absolute function on the sine and cosine. One more rearrangement gives r := reffective * exp( ln( 1 / ( exp(m * ln(abscostheta)) + exp(m * ln(abssintheta)) ) ) / m); which is the form used in the code. In the cases where the sine or cosine are zero (ie on the axes), we must not calculate at all, to avoid log of zero. We simply set r = reffective in those cases. The program has a special feature to speed up the angle of the calculation (theta) so that it moves faster than the angle at which the graph is plotted. With a factor of 3/2, the four corners become 3/2 * 4 = 6 corners, and we obtain a hexagon. examples To produce a nice square, use the parameters: 0.5 First line: lowest value of m. 5.0 Second line: highest value of m. 0.1 Third line: increment in the value of m 1 Fourth line: desired radius of a circle if m = 2. 100 Fifth line: number of steps to take to move around 360 degrees. 1 Sixth line: A factor by which to increase the value of theta. 1 gives a square, 1.5 gives a hexagon. To produce a hexagon transformed into a circle, use the parameters: 1.5 First line: lowest value of m, mlo. 2.0 Second line: highest value of m, mhi. 0.1 Third line: increment in the value of m, mstep. 1.0 Fourth line: desired radius of a circle if m = 2, reffective. 100 Fifth line: number of steps to take to move around 360 degrees. 1.5 Sixth line: A factor by which to increase the value of theta, spinfactor. 1 gives a square, 1.5 gives a hexagon. It is not clear why one has to use the lowest value of n as the same as the theta factor (6th parameter), but it works! (One would have to prove that with these parameters one gets an exact straight hexagon edge.) documentation Inspired by: > Article 7568 in sci.math: > From: pvmg0487@uxa.cso.uiuc.edu > Subject: hexagonal cone function sought > Message-ID: <107700002@uxa.cso.uiuc.edu> > Date: 22 Nov 89 22:25:00 GMT > > I would like to generate a 3-D cone like object, but with a hexagonal > base. Any suggestions as to an appropriate equation? > > Thanks -- Vernon > Article 7578 in sci.math: > From: toms@ncifcrf.gov (Tom Schneider) > Subject: Re: hexagonal cone function sought > Message-ID: <1405@fcs280s.ncifcrf.gov> > Date: 25 Nov 89 01:18:14 GMT > References: <107700002@uxa.cso.uiuc.edu> > Reply-To: toms@fcs260c2.UUCP (Tom Schneider) > Organization: National Cancer Institute, Frederick > Lines: 30 > > In article <107700002@uxa.cso.uiuc.edu> pvmg0487@uxa.cso.uiuc.edu writes: > > > >I would like to generate a 3-D cone like object, but with a hexagonal > >base. Any suggestions as to an appropriate equation? > > > >Thanks -- Vernon > > Well, that's pretty surprising, since just today I was thinking about a > function that does almost exactly what you want! It turns out that the > equation x^n + y^n = r^n is a line (diamond) if n = 1, a circle if n = 2 and > approaches a square as n -> infinity! So all one needs to do is express this > in polar notation, and then scrunch an extra two corners in to get what you > want! > > First, use the form x^n + y^n = rmax^n (to avoid confusion!) and substitute x > = r cos(theta), y = r sin(theta). Divide both sides by r^n, and rearrange to > get r expressed as a function of theta. To get the powers, I had to use a^b > = exp(b*ln(a)). The thing is symmetrical around the 4 quadrants, so I > avoided logs of negative numbers by taking the absolute values of the sine > and cosine functions. Also, at angles of n*pi/2, one gets division by zero, > so just substitute the desired radius. > > I have done this by writing a Pascal program that will do the job. Pretty! > It turns out that to get a hexagon, you have to plot between n=1.5 and n=2 > because of the scrunching. Email me if you want a copy of the program. > > Tom Schneider > National Cancer Institute > Laboratory of Mathematical Biology > Frederick, Maryland 21701-1013 > toms@ncifcrf.gov > From daemon Tue Nov 28 09:34:41 1989 > Return-Path: > Date: Tue, 28 Nov 89 08:35:52 -0600 > From: Paul Vernon McDonald > Message-Id: <8911281435.AA01048@uxa.cso.uiuc.edu> > To: toms@ncifcrf.gov > Subject: Pascal code for hexagon > > Tom, > I'd be most grateful to receive your code, if you are willing to share it. > I curently have a working version of the hexagon, done in piecewise > fashion, but I'd be interested in a generic solution. In fact I plan > to use other shapes in the future, so your code may be of great help. > > Thanks, > > Vernon McDonald > University of Illinois > Department of Kinesiology > Urbana, IL, 61801 > vmcdonald@uiuc.edu > > From toms Tue Nov 28 13:17:17 1989 > To: pvmg0487@uxa.cso.uiuc.edu > Subject: Cisq > > Vernon: > Sure, I wrote the code mostly because of your posting. But actually it > has suddenly become very important to my work (it's a long story...) and so > it is useful to me to have it. I have to brush it up a bit and I will > send it to you. If you have a PostScript printer, then you may also > want the xyplo program, which produces PostScript x-y plotting of data. > This made writing cisq (circle square) easier because I only needed to > create the right numbers and xyplo did the graphics for me. > Tom see also xyplo.p, the Pascal program that produces PostScript x-y plotting graphics. author Thomas Dana Schneider bugs One might also want to produce the hexagon for INCREASING values of n, rather than being confined into the region n=1.5 to 2. It seems that to do this requires that one do a fancy job of warping the square region into the appropriate triangular region. This should be pretty easy with the right afine transformation, but the program doesn't have that feature in it. Fortunately, it is not necessary. technical notes *) (* end module describe.cisq *) (* begin module cisq.const *) wid = 12; (* width of numbers *) dec = 8; (* decimal places of numbers *) firstwid = 3; (* width of first number *) firstdec = 1; (* decimal places of first number *) (* end module cisq.const *) var cisqp: text; (* parameters *) xyin: text; (* input to xyplo *) (* begin module halt *) procedure halt; (* stop the program. the procedure performs a goto to the end of the program. you must have a label: label 1; declared, and also the end of the program must have this label: 1: end. examples are in the module libraries. this is the only goto in the delila system. *) begin writeln(output,' program halt.'); goto 1 end; (* end module halt version = 'delmod 6.57 88 mar 1 tds/gds' *) (* begin module cisq.writeparameters *) procedure writeparameters(var parameters: text; mlo, mhi, mstep, reffective, steps, spinfactor: real); (* write the parameters out *) begin writeln(parameters,'* cisq ',version:4:2); writeln(parameters,'*'); writeln(parameters,'* parameters:'); writeln(parameters,'* mlo = ',mlo:wid:dec); writeln(parameters,'* mhi = ',mhi:wid:dec); writeln(parameters,'* mstep = ',mstep:wid:dec); writeln(parameters,'* reffective = ',reffective:wid:dec); writeln(parameters,'* steps = ',steps:wid:dec); writeln(parameters,'* spinfactor = ',spinfactor:wid:dec); writeln(parameters,'*'); end; (* end module cisq.writeparameters *) (* begin module cisq.readparameters *) procedure readparameters(var parameters: text; var mlo, mhi, mstep, reffective, steps, spinfactor: real); (* read the parameters in, and check them *) procedure nope; (* tell em like it is *) begin (* nope *) writeparameters(output, mlo, mhi, mstep, reffective, steps, spinfactor); halt end; (* nope *) begin reset(parameters); readln(parameters, mlo); readln(parameters, mhi); if mlo > mhi then begin writeln(output,'mlo must be bigger than mhi:'); nope end; readln(parameters, mstep); readln(parameters, reffective); readln(parameters, steps); readln(parameters, spinfactor); if (mlo <= 0) or (mhi <= 0) or (reffective <= 0) or (steps <= 0) or (spinfactor <= 0) or (mhi <= 0) then begin writeln(output,'all parameters must be positive:'); nope end; end; (* end module cisq.readparameters *) (* begin module pic.polrec *) procedure polrec(r,theta: real; var x,y: real); (* convert polar to rectangular coordinates, theta is in radians *) begin x := r*cos(theta); y := r*sin(theta) end; (* end module pic.polrec *) (* begin module cisq.themain *) procedure themain(var cisqp, xyin: text); (* the main procedure of the program *) const pi = 3.1415926535; (* the circumfrence of circle divided by its radius *) var abscostheta: real; (* cosine of theta *) dtheta: real; (* change in theta *) isinteger: boolean; (* true if n is close to an integer by mstep *) m: real; (* power factor in |x|^m + |y|^m = |r|^m *) mhi: real; (* highest value of n *) mlo: real; (* lowest value of n *) mstep: real; (* increments of the value of n *) r: real; (* the radial distance out if n = 2 *) reffective: real; (* the desired final reffective *) sinmstep: real; (* The sine of pi * mstep/2. When sine(pi*m) is close to zero (within sinmstep of zero) then we are near an integer m!! The reason for doing this expensive computation is that the trunc function will not work to indicate where integers are because if one starts with m=0.1 and increments by 0.1, by the time one reaches near 1.0, the value is m=0.9999 so the trunc gives zero. Yet by the time one reaches near 2.0 it is >2.0 on a Sun4. So trunc is inconsistant. Besides, this is elegant! *) abssintheta: real; (* sine of theta *) spinfactor: real; (* factor to make the plotted theta increase by *) steps: real; (* number of steps around the circle *) theta: real; (* angle corresponding to x and y *) twopi: real; (* twice pi *) x: real; (* x coordinate *) y: real; (* y coordinate *) begin rewrite(xyin); readparameters(cisqp, mlo, mhi, mstep, reffective, steps, spinfactor); writeparameters(xyin, mlo, mhi, mstep, reffective, steps, spinfactor); writeparameters(output, mlo, mhi, mstep, reffective, steps, spinfactor); twopi := 2 * pi; dtheta := twopi/steps; m := mlo; (* precalculate the threshold for figuring out when we are near an integer *) sinmstep := sin(pi*mstep/2); (* extra factor of 2 for safety *) while m <= mhi do begin theta := 0.0; if abs(sin(pi*m)) < sinmstep then isinteger := true else isinteger := false; (* xyplo ignores this comment: *) write(xyin,'* m = ',m:wid:dec); if isinteger then writeln(xyin,', an integer') else writeln(xyin); if isinteger or (m = mlo) then begin (* locate the identification tag in the positive quadrant at the intersection of the curve with x=y. This means that 2 |x|^m = reffective^m. Solving for x gives: *) x := reffective / exp(ln(2.0)/m); if m = mlo then write(xyin,'m=',m:firstwid:firstdec) else write(xyin,'m=',round(m):1); writeln(xyin,' ', x:wid:dec, x:wid:dec, (* x = y coordinate *) ' ', 1.0:wid:dec, 1.0:wid:dec); (* dummy coordinates *) end; while theta <= spinfactor*twopi do begin abscostheta := abs(cos(theta)); abssintheta := abs(sin(theta)); (* testing for equality with a real is not so good, but it must be done to work! *) if (abscostheta = 0.0) or (abssintheta = 0.0) then r := reffective else r := reffective * exp( ln( 1 / ( exp(m * ln(abscostheta)) + exp(m * ln(abssintheta)) ) ) / m); polrec(r,theta/spinfactor,x,y); (* slow theta down *) if isinteger then write(xyin,'i') else write(xyin,'r'); writeln(xyin,' ',x:wid:dec,' ',y:wid:dec, ' ',r:wid:dec,' ',theta:wid:dec); theta := theta + dtheta/spinfactor; (* slow dtheta down *) end; m := m + mstep; end; end; (* end module cisq.themain *) begin themain(cisqp, xyin); 1: end.