/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "cisq.p" */ #include /* cisq: circle to square Tom Schneider National Cancer Institute Laboratory of Mathematical Biology Frederick, Maryland 21701-1013 toms@ncifcrf.gov 1989 */ /* end of program */ /* begin module version */ #define 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 */ #define wid 12 /* width of numbers */ #define dec 8 /* decimal places of numbers */ #define firstwid 3 /* width of first number */ #define firstdec 1 /* decimal places of first number */ /* end module cisq.const */ Static _TEXT cisqp; /* parameters */ Static _TEXT xyin; /* input to xyplo */ Static jmp_buf _JL1; /* begin module halt */ Static Void 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. */ printf(" program halt.\n"); longjmp(_JL1, 1); } /* end module halt version = 'delmod 6.57 88 mar 1 tds/gds' */ /* begin module cisq.writeparameters */ Static Void writeparameters(parameters, mlo, mhi, mstep, reffective, steps, spinfactor) _TEXT *parameters; double mlo, mhi, mstep, reffective, steps, spinfactor; { /* write the parameters out */ fprintf(parameters->f, "* cisq %4.2f\n", version); fprintf(parameters->f, "*\n"); fprintf(parameters->f, "* parameters:\n"); fprintf(parameters->f, "* mlo = %*.*f\n", wid, dec, mlo); fprintf(parameters->f, "* mhi = %*.*f\n", wid, dec, mhi); fprintf(parameters->f, "* mstep = %*.*f\n", wid, dec, mstep); fprintf(parameters->f, "* reffective = %*.*f\n", wid, dec, reffective); fprintf(parameters->f, "* steps = %*.*f\n", wid, dec, steps); fprintf(parameters->f, "* spinfactor = %*.*f\n", wid, dec, spinfactor); fprintf(parameters->f, "*\n"); } /* Local variables for readparameters: */ struct LOC_readparameters { double *mlo, *mhi, *mstep, *reffective, *steps, *spinfactor; } ; Local Void nope(LINK) struct LOC_readparameters *LINK; { /* tell em like it is */ _TEXT TEMP; TEMP.f = stdout; *TEMP.name = '\0'; writeparameters(&TEMP, *LINK->mlo, *LINK->mhi, *LINK->mstep, *LINK->reffective, *LINK->steps, *LINK->spinfactor); halt(); } /* nope */ /* end module cisq.writeparameters */ /* begin module cisq.readparameters */ Static Void readparameters(parameters, mlo_, mhi_, mstep_, reffective_, steps_, spinfactor_) _TEXT *parameters; double *mlo_, *mhi_, *mstep_, *reffective_, *steps_, *spinfactor_; { /* read the parameters in, and check them */ struct LOC_readparameters V; V.mlo = mlo_; V.mhi = mhi_; V.mstep = mstep_; V.reffective = reffective_; V.steps = steps_; V.spinfactor = spinfactor_; if (*parameters->name != '\0') { if (parameters->f != NULL) parameters->f = freopen(parameters->name, "r", parameters->f); else parameters->f = fopen(parameters->name, "r"); } else rewind(parameters->f); if (parameters->f == NULL) _EscIO2(FileNotFound, parameters->name); RESETBUF(parameters->f, Char); fscanf(parameters->f, "%lg%*[^\n]", V.mlo); getc(parameters->f); fscanf(parameters->f, "%lg%*[^\n]", V.mhi); getc(parameters->f); if (*V.mlo > *V.mhi) { printf("mlo must be bigger than mhi:\n"); nope(&V); } fscanf(parameters->f, "%lg%*[^\n]", V.mstep); getc(parameters->f); fscanf(parameters->f, "%lg%*[^\n]", V.reffective); getc(parameters->f); fscanf(parameters->f, "%lg%*[^\n]", V.steps); getc(parameters->f); fscanf(parameters->f, "%lg%*[^\n]", V.spinfactor); getc(parameters->f); if (*V.mlo <= 0 || *V.mhi <= 0 || *V.reffective <= 0 || *V.steps <= 0 || *V.spinfactor <= 0 || *V.mhi <= 0) { printf("all parameters must be positive:\n"); nope(&V); } } /* end module cisq.readparameters */ /* begin module pic.polrec */ Static Void polrec(r, theta, x, y) double r, theta, *x, *y; { /* convert polar to rectangular coordinates, theta is in radians */ *x = r * cos(theta); *y = r * sin(theta); } #define pi 3.1415926535 /* the circumfrence of circle divided by its radius */ /* end module pic.polrec */ /* begin module cisq.themain */ Static Void themain(cisqp, xyin) _TEXT *cisqp, *xyin; { /* the main procedure of the program */ double abscostheta; /* cosine of theta */ double dtheta; /* change in theta */ boolean isinteger; /* true if n is close to an integer by mstep */ double m; /* power factor in |x|^m + |y|^m = |r|^m */ double mhi; /* highest value of n */ double mlo; /* lowest value of n */ double mstep; /* increments of the value of n */ double r; /* the radial distance out if n = 2 */ double reffective; /* the desired final reffective */ double sinmstep; /* 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! */ double abssintheta; /* sine of theta */ double spinfactor; /* factor to make the plotted theta increase by */ double steps; /* number of steps around the circle */ double theta; /* angle corresponding to x and y */ double twopi = 2 * pi; /* twice pi */ double x; /* x coordinate */ double y; /* y coordinate */ _TEXT TEMP; if (*xyin->name != '\0') { if (xyin->f != NULL) xyin->f = freopen(xyin->name, "w", xyin->f); else xyin->f = fopen(xyin->name, "w"); } else { if (xyin->f != NULL) rewind(xyin->f); else xyin->f = tmpfile(); } if (xyin->f == NULL) _EscIO2(FileNotFound, xyin->name); SETUPBUF(xyin->f, Char); readparameters(cisqp, &mlo, &mhi, &mstep, &reffective, &steps, &spinfactor); writeparameters(xyin, mlo, mhi, mstep, reffective, steps, spinfactor); TEMP.f = stdout; *TEMP.name = '\0'; writeparameters(&TEMP, mlo, mhi, mstep, reffective, steps, spinfactor); 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) { theta = 0.0; if (fabs(sin(pi * m)) < sinmstep) isinteger = true; else isinteger = false; /* xyplo ignores this comment: */ fprintf(xyin->f, "* m = %*.*f", wid, dec, m); if (isinteger) fprintf(xyin->f, ", an integer\n"); else putc('\n', xyin->f); if (isinteger || m == mlo) { /* 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(log(2.0) / m); if (m == mlo) fprintf(xyin->f, "m=%*.*f", firstwid, firstdec, m); else fprintf(xyin->f, "m=%ld", (long)floor(m + 0.5)); /* x = y coordinate */ fprintf(xyin->f, " %*.*f%*.*f %*.*f%*.*f\n", wid, dec, x, wid, dec, x, wid, dec, 1.0, wid, dec, 1.0); /* dummy coordinates */ } while (theta <= spinfactor * twopi) { abscostheta = fabs(cos(theta)); abssintheta = fabs(sin(theta)); /* testing for equality with a real is not so good, but it must be done to work! */ if (abscostheta == 0.0 || abssintheta == 0.0) r = reffective; else r = reffective * exp( log(1 / (exp(m * log(abscostheta)) + exp(m * log(abssintheta)))) / m); polrec(r, theta / spinfactor, &x, &y); /* slow theta down */ if (isinteger) putc('i', xyin->f); else putc('r', xyin->f); fprintf(xyin->f, " %*.*f %*.*f %*.*f %*.*f\n", wid, dec, x, wid, dec, y, wid, dec, r, wid, dec, theta); theta += dtheta / spinfactor; /* slow dtheta down */ } m += mstep; } } #undef pi /* end module cisq.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; xyin.f = NULL; strcpy(xyin.name, "xyin"); cisqp.f = NULL; strcpy(cisqp.name, "cisqp"); themain(&cisqp, &xyin); _L1: if (cisqp.f != NULL) fclose(cisqp.f); if (xyin.f != NULL) fclose(xyin.f); exit(EXIT_SUCCESS); } /* End. */