This is contour_plot.c in view mode; [Download] [Up]
/* contour_plot.f -- translated by f2c (version of 22 January 1991 19:25:10).
You must link the resulting object file with the libraries:
-lF77 -lI77 -lm -lc (in that order)
*/
#import <appkit/color.h>
#import "f2c.h"
#import <dpsclient/wraps.h>
/* Common Block Declarations */
struct {
real xt[1000], yt[1000], zt[1000];
integer ia[3000];
} temp1_;
#define temp1_1 temp1_
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int contour_(jdim, kdim, nj, nk, x, y, f, ncont, acont,
ppxunit, ppyunit, colorMap)
integer *jdim, *kdim, *nj, *nk, *colorMap;
real *x, *y, *f;
integer *ncont;
real *acont;
float *ppxunit, *ppyunit;
{
/* System generated locals */
integer x_dim1, x_offset, y_dim1, y_offset, f_dim1, f_offset;
/* Local variables */
extern /* Subroutine */ int con2l_();
/* Parameter adjustments */
--acont;
f_dim1 = *jdim;
f_offset = f_dim1 + 1;
f -= f_offset;
y_dim1 = *jdim;
y_offset = y_dim1 + 1;
y -= y_offset;
x_dim1 = *jdim;
x_offset = x_dim1 + 1;
x -= x_offset;
/* Function Body */
con2l_(jdim, kdim, &c__1, nj, &c__1, nk, &x[x_offset], &y[y_offset], &f[
f_offset], ncont, &acont[1], ppxunit, ppyunit, colorMap);
return 0;
} /* contour_ */
/* Subroutine */ int cline2_(cont, np, x, y,
ppxunit, ppyunit, colorMap, icont, ncont)
real *cont;
integer *np, *colorMap;
int icont, *ncont;
real *x, *y;
float *ppxunit, *ppyunit;
{
int i;
float pattern0[] = {}; /* solid */
float pattern1[] = {4.0, 4.0}; /* dash */
NXColor color;
float r,g,b;
float h,s,br;
int one_third;
/* Parameter adjustments */
--y;
--x;
/* PSsetlinewidth(0.0); set outside ! */
/* PSsetgray(0); set outside ! */
switch(*colorMap) {
case 0:
PSsetdash(pattern0, 0, 0.0);
break;
case 1:
if(*cont >= 0.0)PSsetdash(pattern0, 0, 0.0);
if(*cont < 0.0) PSsetdash(pattern1, 2, 0.0);
break;
case 2:
PSsetdash(pattern0, 0, 0.0);
one_third = *ncont/3;
r = 0;
g = 0;
b = 0;
h = 0;
s = 0;
br = 0;
/* if(icont <= one_third)
{
b = 1.-(float)icont/one_third;
g = (float)icont/one_third;
r = b*g;
}
if(icont > one_third && icont <= 2*one_third)
{
b = ((float)icont - (float)one_third)/one_third;
g = 1. - ((float)icont - (float)one_third)/one_third;
r = g*b;
}
if(icont > 2*one_third && icont <= *ncont)
{
r = 1.-((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
b = ((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
g = r*b;
}
color = NXConvertRGBToColor(r,g,b);
*/
h = 0.6666*icont/(float)*ncont;
s = 1.0;
br = 1.0;
color = NXConvertHSBToColor(h,s,br);
NXSetColor(color);
break;
case 3:
if(*cont >= 0.0)PSsetdash(pattern0, 0, 0.0);
if(*cont < 0.0) PSsetdash(pattern1, 2, 0.0);
one_third = *ncont/3;
r = 0;
g = 0;
b = 0;
h = 0;
s = 0;
br = 0;
/* if(icont <= one_third)
{
b = 1.-(float)icont/one_third;
g = (float)icont/one_third;
r = b*g;
}
if(icont > one_third && icont <= 2*one_third)
{
b = ((float)icont - (float)one_third)/one_third;
g = 1. - ((float)icont - (float)one_third)/one_third;
r = g*b;
}
if(icont > 2*one_third && icont <= *ncont)
{
r = 1.-((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
b = ((float)icont-2.*(float)one_third)/((float)*ncont-2.*one_third);
g = r*b;
}
color = NXConvertRGBToColor(r,g,b);
*/
h = 0.6666*icont/(float)*ncont;
s = 1.0;
br = 1.0;
color = NXConvertHSBToColor(h,s,br);
NXSetColor(color);
break;
}
PSnewpath();
PSmoveto(x[1]*(*ppxunit), y[1]*(*ppyunit));
for (i=2; i <= *np; i++){
PSlineto(x[i]*(*ppxunit), y[i]*(*ppyunit));
}
PSstroke();
/* Function Body */
return 0;
} /* cline2_ */
/* ----------------------------------------------------------------------- */
/* Subroutine */ int con2l_(idim, jdim, is, ie, js, je, x, y, f, ncont, acont,
ppxunit, ppyunit, colorMap)
integer *idim, *jdim, *is, *ie, *js, *je, *colorMap;
real *x, *y, *f;
integer *ncont;
real *acont;
float *ppxunit, *ppyunit;
{
/* System generated locals */
integer x_dim1, x_offset, y_dim1, y_offset, f_dim1, f_offset, i__1, i__2,
i__3, i__4;
/* Local variables */
static integer ibeg, jbeg, idir;
static real cont;
static integer i, j, iaend, icont;
extern /* Subroutine */ int cline2_();
static integer ij, np;
static real wx, wy;
static integer ignext, lnstrt, iae, iia, iee, jee, imb, ima;
static real fij;
static integer iss, jss;
static real xyf, wxx, wyy;
/* Calculate contour lines for the function F in the region IS to IE, JS
to*/
/* JE. X,Y coordinates corresponding to the grid points are in arrays X
and*/
/* Y. */
/* Common block TEMP1 contains scratch arrays XT, YT, ZT, and IA, and may
be*/
/* used elsewhere. This subroutine is supposed to be able to recover fro
m*/
/* filling either array. */
/* One little check. If IS=IE or JS=JE, return with no contour lines.
*/
/* Parameter adjustments */
--acont;
f_dim1 = *idim;
f_offset = f_dim1 + 1;
f -= f_offset;
y_dim1 = *idim;
y_offset = y_dim1 + 1;
y -= y_offset;
x_dim1 = *idim;
x_offset = x_dim1 + 1;
x -= x_offset;
/* Function Body */
if (*is == *ie || *js == *je) {
return 0;
}
/* Zero the marker array. */
for (iia = 1; iia <= 3000; ++iia) {
temp1_1.ia[iia - 1] = 0;
/* L100: */
}
/* LNSTRT=1 means we're starting a new line. */
lnstrt = 1;
ignext = 0;
/* Loop through each contour level. */
i__1 = *ncont;
for (icont = 1; icont <= i__1; ++icont) {
cont = acont[icont];
/* The IS-IE,JS-JE region may have to be subdivided because of IA spa
ce. This*/
/* may have a detrimental effect on labelling of the contour lines.
Move*/
/* IS,IE,JS,JE to new variables. */
iss = *is;
iee = *ie;
jss = *js;
jee = *je;
/* Flag points in IA where the the function increases through the con
tour*/
/* level, not including the boundaries. This is so we have a list of
at least*/
/* one point on each contour line that doesn't intersect a boundary.
*/
L110:
iae = 0;
i__2 = jee - 1;
for (j = jss + 1; j <= i__2; ++j) {
imb = 0;
iaend = iae;
i__3 = iee;
for (i = iss; i <= i__3; ++i) {
if (f[i + j * f_dim1] <= cont) {
imb = 1;
} else if (imb == 1) {
++iae;
temp1_1.ia[iae - 1] = i * 1000 + j;
imb = 0;
/* Check if the IA array is full. If so, the subdividing
algorithm goes like*/
/* this: if we've marked at least one J row, drop back t
o the last completed*/
/* J and call that the region. If we haven't even finish
ed one J row, our*/
/* region just extends to this I location. */
if (iae == 3000) {
if (j > jss + 1) {
iae = iaend;
jee = j;
} else {
/* Computing MIN */
i__4 = j + 1;
jee = min(i__4,jee);
iee = i;
}
goto L210;
}
}
/* L200: */
}
}
/* Search along the boundaries for contour line starts. IMA tells w
hich */
/* boundary of the region we're on. */
L210:
ima = 1;
imb = 0;
ibeg = iss - 1;
jbeg = jss;
L1:
++ibeg;
if (ibeg == iee) {
ima = 2;
}
goto L5;
L2:
++jbeg;
if (jbeg == jee) {
ima = 3;
}
goto L5;
L3:
--ibeg;
if (ibeg == iss) {
ima = 4;
}
goto L5;
L4:
--jbeg;
if (jbeg == jss) {
ima = 5;
}
L5:
if (f[ibeg + jbeg * f_dim1] > cont) {
goto L7;
}
imb = 1;
L6:
switch ((int)ima) {
case 1: goto L1;
case 2: goto L2;
case 3: goto L3;
case 4: goto L4;
case 5: goto L91;
}
L7:
if (imb != 1) {
goto L6;
}
/* Got a start point. */
imb = 0;
i = ibeg;
j = jbeg;
fij = f[ibeg + jbeg * f_dim1];
/* Round the corner if necessary. */
switch ((int)ima) {
case 1: goto L21;
case 2: goto L11;
case 3: goto L12;
case 4: goto L13;
case 5: goto L51;
}
L11:
if (j != jss) {
goto L31;
}
goto L21;
L12:
if (i != iee) {
goto L41;
}
goto L31;
L13:
if (j != jee) {
goto L51;
}
goto L41;
/* This is how we start in the middle of the region, using IA. */
L10:
i = temp1_1.ia[iia - 1] / 1000;
j = temp1_1.ia[iia - 1] - i * 1000;
fij = f[i + j * f_dim1];
temp1_1.ia[iia - 1] = 0;
goto L21;
/*
4 */
/* Look different directions to see which way the contour line went:
1-|-3*/
/*
2 */
L20:
++j;
L21:
--i;
if (i < iss) {
goto L90;
}
idir = 1;
if (f[i + j * f_dim1] <= cont) {
goto L52;
}
fij = f[i + j * f_dim1];
goto L31;
L30:
--i;
L31:
--j;
if (j < jss) {
goto L90;
}
idir = 2;
if (f[i + j * f_dim1] <= cont) {
goto L60;
}
fij = f[i + j * f_dim1];
goto L41;
L40:
--j;
L41:
++i;
if (i > iee) {
goto L90;
}
idir = 3;
if (f[i + j * f_dim1] <= cont) {
goto L60;
}
fij = f[i + j * f_dim1];
goto L51;
L50:
++i;
L51:
++j;
idir = 4;
if (j > jee) {
goto L90;
}
if (f[i + j * f_dim1] <= cont) {
goto L60;
}
fij = f[i + j * f_dim1];
goto L21;
/* Wipe this point out of IA if it's in the list. */
L52:
if (iae == 0) {
goto L60;
}
ij = i * 1000 + j + 1000;
i__3 = iae;
for (iia = 1; iia <= i__3; ++iia) {
if (temp1_1.ia[iia - 1] == ij) {
temp1_1.ia[iia - 1] = 0;
goto L60;
}
/* L300: */
}
/* Do interpolation for X,Y coordinates. */
L60:
xyf = (cont - f[i + j * f_dim1]) / (fij - f[i + j * f_dim1]);
/* This tests for a contour point coinciding with a grid point. In t
his case*/
/* the contour routine comes up with the same physical coordinate twi
ce. If*/
/* If we don't trap it, it can (in some cases significantly) increase
the*/
/* number of points in a contour line. Also, if this happens on the
first*/
/* point in a line, the second point could be misinterpreted as the e
nd of a*/
/* (circling) contour line. */
if (xyf == (float)0.) {
++ignext;
}
switch ((int)idir) {
case 1: goto L61;
case 2: goto L62;
case 3: goto L63;
case 4: goto L64;
}
L61:
wxx = x[i + j * x_dim1] + xyf * (x[i + 1 + j * x_dim1] - x[i + j *
x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i + 1 + j * y_dim1] - y[i + j *
y_dim1]);
goto L65;
L62:
wxx = x[i + j * x_dim1] + xyf * (x[i + (j + 1) * x_dim1] - x[i + j *
x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i + (j + 1) * y_dim1] - y[i + j *
y_dim1]);
goto L65;
L63:
wxx = x[i + j * x_dim1] + xyf * (x[i - 1 + j * x_dim1] - x[i + j *
x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i - 1 + j * y_dim1] - y[i + j *
y_dim1]);
goto L65;
L64:
wxx = x[i + j * x_dim1] + xyf * (x[i + (j - 1) * x_dim1] - x[i + j *
x_dim1]);
wyy = y[i + j * y_dim1] + xyf * (y[i + (j - 1) * y_dim1] - y[i + j *
y_dim1]);
/* Figure out what to do with this point. */
L65:
if (lnstrt != 1) {
goto L66;
}
/* This is the first point in a contour line. */
np = 1;
temp1_1.xt[np - 1] = wxx;
temp1_1.yt[np - 1] = wyy;
wx = wxx;
wy = wyy;
lnstrt = 0;
goto L67;
/* Add a point to this line. Check for duplicate point first. */
L66:
if (ignext == 2) {
if (wxx == temp1_1.xt[np - 1] && wyy == temp1_1.yt[np - 1]) {
ignext = 0;
goto L67;
} else {
ignext = 1;
}
}
++np;
temp1_1.xt[np - 1] = wxx;
temp1_1.yt[np - 1] = wyy;
/* See if the temporary array xT is full. */
if (np == 1000) {
cline2_(&cont, &np, temp1_1.xt, temp1_1.yt,
ppxunit, ppyunit, colorMap, icont, ncont);
temp1_1.xt[0] = temp1_1.xt[np - 1];
temp1_1.yt[0] = temp1_1.yt[np - 1];
np = 1;
}
/* Check to see if we're back to the initial point. */
if (wxx != wx) {
goto L67;
}
if (wyy == wy) {
goto L90;
}
/* Nope. Search for the next point on this line. */
L67:
switch ((int)idir) {
case 1: goto L50;
case 2: goto L20;
case 3: goto L30;
case 4: goto L40;
}
/* This is the end of a contour line. After this we'll start a new l
ine.*/
L90:
lnstrt = 1;
ignext = 0;
cline2_(&cont, &np, temp1_1.xt, temp1_1.yt,
ppxunit, ppyunit, colorMap, icont, ncont);
/* If we're not done looking along the boundaries, go look there some
more.*/
if (ima != 5) {
goto L6;
}
/* Otherwise, get the next start out of IA. */
L91:
if (iae != 0) {
i__3 = iae;
for (iia = 1; iia <= i__3; ++iia) {
if (temp1_1.ia[iia - 1] != 0) {
goto L10;
}
/* L400: */
}
}
/* And if there are no more of these, we're done with this region. I
f we've*/
/* subdivided, update the region pointers and go back for more. */
if (iee == *ie) {
if (jee != *je) {
jss = jee;
jee = *je;
goto L110;
}
} else {
iss = iee;
iee = *ie;
goto L110;
}
/* Loop back for the next contour level. */
/* L1000: */
}
return 0;
} /* con2l_ */
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.