00001 /* srface.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Common Block Declarations */
00009 
00010 struct srfblk_1_ {
00011     integer limu[1024], liml[1024];
00012     real cl[41];
00013     integer ncl, ll;
00014     real fact;
00015     integer irot, ndrz, nupper, nrswt;
00016     real bigd, umin, umax, vmin, vmax, rzero;
00017     integer ioffp, nspval;
00018     real spval, bigest;
00019 };
00020 
00021 #define srfblk_1 (*(struct srfblk_1_ *) &srfblk_)
00022 
00023 struct {
00024     real xxmin, xxmax, yymin, yymax, zzmin, zzmax, delcrt, eyex, eyey, eyez;
00025 } pwrz1s_;
00026 
00027 #define pwrz1s_1 pwrz1s_
00028 
00029 struct srfip1_1_ {
00030     integer ifr, istp, irots, idrx, idry, idrz, iupper, iskirt, ncla;
00031     real theta, hskirt, chi, clo, cinc;
00032     integer ispval;
00033 };
00034 
00035 #define srfip1_1 (*(struct srfip1_1_ *) &srfip1_)
00036 
00037 /* Initialized data */
00038 
00039 struct {
00040     integer fill_1[2095];
00041     integer e_2;
00042     integer fill_3[6];
00043     integer e_4;
00044     integer fill_5[1];
00045     real e_6;
00046     integer fill_7[1];
00047     } srfblk_ = { {0}, 0, {0}, 0, {0}, 0.f };
00048 
00049 struct {
00050     integer e_1[9];
00051     real e_2[5];
00052     integer e_3;
00053     } srfip1_ = { 1, 0, 0, 1, 1, 0, 0, 0, 6, .02f, 0.f, 0.f, 0.f, 0.f, -999 };
00054 
00055 
00056 /* Table of constant values */
00057 
00058 static real c_b2 = 0.f;
00059 static real c_b3 = 1.f;
00060 static real c_b7 = 1024.f;
00061 static integer c__1 = 1;
00062 static integer c__40 = 40;
00063 static integer c__0 = 0;
00064 static real c_b128 = 10.f;
00065 static integer c__2 = 2;
00066 
00067 /* ======================================================================= */
00068 
00069 
00070 
00071 /* Subroutine */ int srface_(real *x, real *y, real *z__, integer *m, integer 
00072         *mx, integer *nx, integer *ny, real *s, real *stereo)
00073 {
00074     /* Initialized data */
00075 
00076     static integer jf = 1;
00077     static integer if__ = 1;
00078     static integer ly = 2;
00079     static integer lx = 2;
00080     static integer icnst = 0;
00081 
00082     /* System generated locals */
00083     integer z_dim1, z_offset, m_dim2, m_offset, i__1, i__2, i__3, i__4, i__5, 
00084             i__6, i__7, i__8, i__9, i__10;
00085 
00086     /* Local variables */
00087     static integer ipic, npic, ipli, jplj;
00088     static real ster, poix, poiy, poiz, xeye;
00089     static integer mmxx, nnxx;
00090     static real yeye;
00091     static integer nnyy;
00092     static real zeye, ynow, xnow, sign1;
00093     static integer i__, j, k, l;
00094     extern /* Subroutine */ int frame_(void);
00095     static real hight;
00096     extern /* Subroutine */ int clset_(real *, integer *, integer *, integer *
00097             , real *, real *, real *, integer *, integer *, real *, integer *,
00098              integer *, integer *, real *, real *);
00099     static real width;
00100     extern /* Subroutine */ int draws_(integer *, integer *, integer *, 
00101             integer *, integer *, integer *);
00102     static integer jpass, ipass;
00103     static real d1, d2;
00104     extern /* Subroutine */ int trn32s_(real *, real *, real *, real *, real *
00105             , real *, integer *);
00106     static real dummy;
00107     static integer nxstp, nystp, ii, jj, li, mi, in, jn, ni, lj;
00108     static real dx, dy;
00109     static integer mj, nj;
00110     extern /* Subroutine */ int srfabd_(void);
00111     static real ctheta, rx, ry, rz, ut, vt, qu, qv, ru, zz, rv;
00112     extern /* Subroutine */ int ctcell_(real *, integer *, integer *, integer 
00113             *, integer *, integer *, integer *);
00114     static real stheta;
00115     static integer nxpass, nypass;
00116     static real ux1, vx1, ux2, vx2, uy1, vy1, uy2, vy2, dif, agl;
00117     static integer nla, mxf[2], myf[2];
00118     extern /* Subroutine */ int set_(real *, real *, real *, real *, real *, 
00119             real *, real *, real *, integer *);
00120     static integer mxj[2], myj[2], mxs[2], mys[2], nxp1, nyp1;
00121 
00122 
00123 /*  Surface plotting package from NCAR -- the only high level NCAR */
00124 /*  routine in this library at present (Aug 17, 1990). */
00125 
00126 /*cc      DIMENSION       X(NX)      ,Y(NY)      ,Z(MX,NY)   ,M(2,NX,NY) ,
00127 */
00128 /* cc     1                S(6) */
00129     /* Parameter adjustments */
00130     --x;
00131     m_dim2 = *nx;
00132     m_offset = (m_dim2 + 1 << 1) + 1;
00133     m -= m_offset;
00134     z_dim1 = *mx;
00135     z_offset = z_dim1 + 1;
00136     z__ -= z_offset;
00137     --y;
00138     --s;
00139 
00140     /* Function Body */
00141 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00142  */
00143     srfabd_();
00144     set_(&c_b2, &c_b3, &c_b2, &c_b3, &c_b3, &c_b7, &c_b3, &c_b7, &c__1);
00145     srfblk_1.bigest = 1e38f;
00146 /* CC      BIGEST = R1MACH(2) */
00147     mmxx = *mx;
00148     nnxx = *nx;
00149     nnyy = *ny;
00150     ster = *stereo;
00151     nxp1 = nnxx + 1;
00152     nyp1 = nnyy + 1;
00153     nla = srfip1_1.ncla;
00154     srfblk_1.nspval = srfip1_1.ispval;
00155     srfblk_1.ndrz = srfip1_1.idrz;
00156     if (srfip1_1.idrz != 0) {
00157         clset_(&z__[z_offset], &mmxx, &nnxx, &nnyy, &srfip1_1.chi, &
00158                 srfip1_1.clo, &srfip1_1.cinc, &nla, &c__40, srfblk_1.cl, &
00159                 srfblk_1.ncl, &icnst, &srfblk_1.ioffp, &srfblk_1.spval, &
00160                 srfblk_1.bigest);
00161     }
00162     if (srfip1_1.idrz != 0) {
00163         srfblk_1.ndrz = 1 - icnst;
00164     }
00165     stheta = sin(ster * srfip1_1.theta);
00166     ctheta = cos(ster * srfip1_1.theta);
00167     rx = s[1] - s[4];
00168     ry = s[2] - s[5];
00169     rz = s[3] - s[6];
00170     d1 = sqrt(rx * rx + ry * ry + rz * rz);
00171     d2 = sqrt(rx * rx + ry * ry);
00172     dx = 0.f;
00173     dy = 0.f;
00174     if (*stereo == 0.f) {
00175         goto L20;
00176     }
00177     d1 = d1 * *stereo * srfip1_1.theta;
00178     if (d2 > 0.f) {
00179         goto L10;
00180     }
00181     dx = d1;
00182     goto L20;
00183 L10:
00184     agl = atan2(rx, -ry);
00185     dx = d1 * cos(agl);
00186     dy = d1 * sin(agl);
00187 L20:
00188     srfblk_1.irot = srfip1_1.irots;
00189     npic = 1;
00190     if (ster != 0.f) {
00191         npic = 2;
00192     }
00193     srfblk_1.fact = 1.f;
00194     if (srfblk_1.nrswt != 0) {
00195         srfblk_1.fact = srfblk_1.rzero / d1;
00196     }
00197     if (srfip1_1.istp == 0 && ster != 0.f) {
00198         srfblk_1.irot = 1;
00199     }
00200     i__1 = npic;
00201     for (ipic = 1; ipic <= i__1; ++ipic) {
00202         srfblk_1.nupper = srfip1_1.iupper;
00203         if (srfip1_1.ifr < 0) {
00204             frame_();
00205         }
00206 
00207 /* SET UP MAPING FROM FLOATING POINT 3-SPACE TO CRT SPACE. */
00208 
00209         sign1 = (real) ((ipic << 1) - 3);
00210         pwrz1s_1.eyex = s[1] + sign1 * dx;
00211         poix = s[4] + sign1 * dx;
00212         pwrz1s_1.eyey = s[2] + sign1 * dy;
00213         poiy = s[5] + sign1 * dy;
00214         pwrz1s_1.eyez = s[3];
00215         poiz = s[6];
00216         srfblk_1.ll = 0;
00217         xeye = pwrz1s_1.eyex;
00218         yeye = pwrz1s_1.eyey;
00219         zeye = pwrz1s_1.eyez;
00220         trn32s_(&poix, &poiy, &poiz, &xeye, &yeye, &zeye, &c__0);
00221         srfblk_1.ll = ipic + (srfip1_1.istp << 1) + 3;
00222         if (ster == 0.f) {
00223             srfblk_1.ll = 1;
00224         }
00225         if (srfblk_1.nrswt != 0) {
00226             goto L100;
00227         }
00228         pwrz1s_1.xxmin = x[1];
00229         pwrz1s_1.xxmax = x[nnxx];
00230         pwrz1s_1.yymin = y[1];
00231         pwrz1s_1.yymax = y[nnyy];
00232         srfblk_1.umin = srfblk_1.bigest;
00233         srfblk_1.vmin = srfblk_1.bigest;
00234         pwrz1s_1.zzmin = srfblk_1.bigest;
00235         srfblk_1.umax = -srfblk_1.umin;
00236         srfblk_1.vmax = -srfblk_1.vmin;
00237         pwrz1s_1.zzmax = -pwrz1s_1.zzmin;
00238         i__2 = nnyy;
00239         for (j = 1; j <= i__2; ++j) {
00240             i__3 = nnxx;
00241             for (i__ = 1; i__ <= i__3; ++i__) {
00242                 zz = z__[i__ + j * z_dim1];
00243                 if (srfblk_1.ioffp == 1 && zz == srfblk_1.spval) {
00244                     goto L30;
00245                 }
00246                 pwrz1s_1.zzmax = dmax(pwrz1s_1.zzmax,zz);
00247                 pwrz1s_1.zzmin = dmin(pwrz1s_1.zzmin,zz);
00248                 trn32s_(&x[i__], &y[j], &z__[i__ + j * z_dim1], &ut, &vt, &
00249                         dummy, &c__1);
00250                 srfblk_1.umax = dmax(srfblk_1.umax,ut);
00251                 srfblk_1.umin = dmin(srfblk_1.umin,ut);
00252                 srfblk_1.vmax = dmax(srfblk_1.vmax,vt);
00253                 srfblk_1.vmin = dmin(srfblk_1.vmin,vt);
00254 L30:
00255                 ;
00256             }
00257 /* L40: */
00258         }
00259         if (srfip1_1.iskirt != 1) {
00260             goto L70;
00261         }
00262         nxstp = nnxx - 1;
00263         nystp = nnyy - 1;
00264         i__2 = nnyy;
00265         i__3 = nystp;
00266         for (j = 1; i__3 < 0 ? j >= i__2 : j <= i__2; j += i__3) {
00267             i__4 = nnxx;
00268             i__5 = nxstp;
00269             for (i__ = 1; i__5 < 0 ? i__ >= i__4 : i__ <= i__4; i__ += i__5) {
00270                 trn32s_(&x[i__], &y[j], &srfip1_1.hskirt, &ut, &vt, &dummy, &
00271                         c__1);
00272                 srfblk_1.umax = dmax(srfblk_1.umax,ut);
00273                 srfblk_1.umin = dmin(srfblk_1.umin,ut);
00274                 srfblk_1.vmax = dmax(srfblk_1.vmax,vt);
00275                 srfblk_1.vmin = dmin(srfblk_1.vmin,vt);
00276 /* L50: */
00277             }
00278 /* L60: */
00279         }
00280 L70:
00281         width = srfblk_1.umax - srfblk_1.umin;
00282         hight = srfblk_1.vmax - srfblk_1.vmin;
00283         dif = (width - hight) * .5f;
00284         if (dif < 0.f) {
00285             goto L80;
00286         } else if (dif == 0) {
00287             goto L100;
00288         } else {
00289             goto L90;
00290         }
00291 L80:
00292         srfblk_1.umin += dif;
00293         srfblk_1.umax -= dif;
00294         goto L100;
00295 L90:
00296         srfblk_1.vmin -= dif;
00297         srfblk_1.vmax += dif;
00298 L100:
00299         xeye = pwrz1s_1.eyex;
00300         yeye = pwrz1s_1.eyey;
00301         zeye = pwrz1s_1.eyez;
00302         trn32s_(&poix, &poiy, &poiz, &xeye, &yeye, &zeye, &c__0);
00303         i__3 = nnyy;
00304         for (j = 1; j <= i__3; ++j) {
00305             i__2 = nnxx;
00306             for (i__ = 1; i__ <= i__2; ++i__) {
00307                 trn32s_(&x[i__], &y[j], &z__[i__ + j * z_dim1], &ut, &vt, &
00308                         dummy, &c__1);
00309                 m[(i__ + j * m_dim2 << 1) + 1] = ut;
00310                 m[(i__ + j * m_dim2 << 1) + 2] = vt;
00311 /* L110: */
00312             }
00313 /* L120: */
00314         }
00315 
00316 /* INITIALIZE UPPER AND LOWER VISIBILITY ARRAYS */
00317 
00318         for (k = 1; k <= 1024; ++k) {
00319             srfblk_1.limu[k - 1] = 0;
00320             srfblk_1.liml[k - 1] = 1024;
00321 /* L130: */
00322         }
00323 
00324 /* FIND ORDER TO DRAW LINES */
00325 
00326         nxpass = 1;
00327         if (s[1] >= x[nnxx]) {
00328             goto L160;
00329         }
00330         if (s[1] <= x[1]) {
00331             goto L170;
00332         }
00333         i__3 = nnxx;
00334         for (i__ = 2; i__ <= i__3; ++i__) {
00335             lx = i__;
00336             if (s[1] <= x[i__]) {
00337                 goto L150;
00338             }
00339 /* L140: */
00340         }
00341 L150:
00342         mxs[0] = lx - 1;
00343         mxj[0] = -1;
00344         mxf[0] = 1;
00345         mxs[1] = lx;
00346         mxj[1] = 1;
00347         mxf[1] = nnxx;
00348         nxpass = 2;
00349         goto L180;
00350 L160:
00351         mxs[0] = nnxx;
00352         mxj[0] = -1;
00353         mxf[0] = 1;
00354         goto L180;
00355 L170:
00356         mxs[0] = 1;
00357         mxj[0] = 1;
00358         mxf[0] = nnxx;
00359 L180:
00360         nypass = 1;
00361         if (s[2] >= y[nnyy]) {
00362             goto L210;
00363         }
00364         if (s[2] <= y[1]) {
00365             goto L220;
00366         }
00367         i__3 = nnyy;
00368         for (j = 2; j <= i__3; ++j) {
00369             ly = j;
00370             if (s[2] <= y[j]) {
00371                 goto L200;
00372             }
00373 /* L190: */
00374         }
00375 L200:
00376         mys[0] = ly - 1;
00377         myj[0] = -1;
00378         myf[0] = 1;
00379         mys[1] = ly;
00380         myj[1] = 1;
00381         myf[1] = nnyy;
00382         nypass = 2;
00383         goto L230;
00384 L210:
00385         mys[0] = nnyy;
00386         myj[0] = -1;
00387         myf[0] = 1;
00388         goto L230;
00389 L220:
00390         mys[0] = 1;
00391         myj[0] = 1;
00392         myf[0] = nnyy;
00393 
00394 /* PUT ON SKIRT ON FRONT SIDE IF WANTED */
00395 
00396 L230:
00397         if (nxpass == 2 && nypass == 2) {
00398             goto L490;
00399         }
00400         if (srfip1_1.iskirt == 0) {
00401             goto L290;
00402         }
00403         in = mxs[0];
00404         if__ = mxf[0];
00405         jn = mys[0];
00406         jf = myf[0];
00407         if (nypass != 1) {
00408             goto L260;
00409         }
00410         trn32s_(&x[1], &y[jn], &srfip1_1.hskirt, &ux1, &vx1, &dummy, &c__1);
00411         trn32s_(&x[nnxx], &y[jn], &srfip1_1.hskirt, &ux2, &vx2, &dummy, &c__1)
00412                 ;
00413         qu = (ux2 - ux1) / (x[nnxx] - x[1]);
00414         qv = (vx2 - vx1) / (x[nnxx] - x[1]);
00415         ynow = y[jn];
00416         i__3 = nnxx;
00417         for (i__ = 1; i__ <= i__3; ++i__) {
00418             trn32s_(&x[i__], &ynow, &srfip1_1.hskirt, &ru, &rv, &dummy, &c__1)
00419                     ;
00420             i__2 = (integer) ru;
00421             i__5 = (integer) rv;
00422             draws_(&i__2, &i__5, &m[(i__ + jn * m_dim2 << 1) + 1], &m[(i__ + 
00423                     jn * m_dim2 << 1) + 2], &c__1, &c__0);
00424 /* L240: */
00425         }
00426         i__3 = (integer) ux1;
00427         i__2 = (integer) vx1;
00428         i__5 = (integer) ux2;
00429         i__4 = (integer) vx2;
00430         draws_(&i__3, &i__2, &i__5, &i__4, &c__1, &c__1);
00431         if (srfip1_1.idry != 0) {
00432             goto L260;
00433         }
00434         i__3 = nnxx;
00435         for (i__ = 2; i__ <= i__3; ++i__) {
00436             draws_(&m[(i__ - 1 + jn * m_dim2 << 1) + 1], &m[(i__ - 1 + jn * 
00437                     m_dim2 << 1) + 2], &m[(i__ + jn * m_dim2 << 1) + 1], &m[(
00438                     i__ + jn * m_dim2 << 1) + 2], &c__1, &c__1);
00439 /* L250: */
00440         }
00441 L260:
00442         if (nxpass != 1) {
00443             goto L290;
00444         }
00445         trn32s_(&x[in], &y[1], &srfip1_1.hskirt, &uy1, &vy1, &dummy, &c__1);
00446         trn32s_(&x[in], &y[nnyy], &srfip1_1.hskirt, &uy2, &vy2, &dummy, &c__1)
00447                 ;
00448         qu = (uy2 - uy1) / (y[nnyy] - y[1]);
00449         qv = (vy2 - vy1) / (y[nnyy] - y[1]);
00450         xnow = x[in];
00451         i__3 = nnyy;
00452         for (j = 1; j <= i__3; ++j) {
00453             trn32s_(&xnow, &y[j], &srfip1_1.hskirt, &ru, &rv, &dummy, &c__1);
00454             i__2 = (integer) ru;
00455             i__5 = (integer) rv;
00456             draws_(&i__2, &i__5, &m[(in + j * m_dim2 << 1) + 1], &m[(in + j * 
00457                     m_dim2 << 1) + 2], &c__1, &c__0);
00458 /* L270: */
00459         }
00460         i__3 = (integer) uy1;
00461         i__2 = (integer) vy1;
00462         i__5 = (integer) uy2;
00463         i__4 = (integer) vy2;
00464         draws_(&i__3, &i__2, &i__5, &i__4, &c__1, &c__1);
00465         if (srfip1_1.idrx != 0) {
00466             goto L290;
00467         }
00468         i__3 = nnyy;
00469         for (j = 2; j <= i__3; ++j) {
00470             draws_(&m[(in + (j - 1) * m_dim2 << 1) + 1], &m[(in + (j - 1) * 
00471                     m_dim2 << 1) + 2], &m[(in + j * m_dim2 << 1) + 1], &m[(in 
00472                     + j * m_dim2 << 1) + 2], &c__1, &c__1);
00473 /* L280: */
00474         }
00475 
00476 /* PICK PROPER ALGORITHM */
00477 
00478 L290:
00479         li = mxj[0];
00480         mi = mxs[0] - li;
00481         ni = (i__3 = mi - mxf[0], abs(i__3));
00482         lj = myj[0];
00483         mj = mys[0] - lj;
00484         nj = (i__3 = mj - myf[0], abs(i__3));
00485 
00486 /* WHEN LINE OF SIGHT IS NEARER TO PARALLEL TO THE X AXIS, */
00487 /* HAVE J LOOP OUTER-MOST, OTHERWISE HAVE I LOOP OUTER-MOST. */
00488 
00489         if (dabs(rx) <= dabs(ry)) {
00490             goto L360;
00491         }
00492         if (srfip1_1.iskirt != 0 || nypass != 1) {
00493             goto L310;
00494         }
00495         i__ = mxs[0];
00496         i__3 = nnyy;
00497         for (j = 2; j <= i__3; ++j) {
00498             draws_(&m[(i__ + (j - 1) * m_dim2 << 1) + 1], &m[(i__ + (j - 1) * 
00499                     m_dim2 << 1) + 2], &m[(i__ + j * m_dim2 << 1) + 1], &m[(
00500                     i__ + j * m_dim2 << 1) + 2], &c__0, &c__1);
00501 /* L300: */
00502         }
00503 L310:
00504         i__3 = nnxx;
00505         for (ii = 1; ii <= i__3; ++ii) {
00506             i__ = mi + ii * li;
00507             ipli = i__ + li;
00508             if (nypass == 1) {
00509                 goto L320;
00510             }
00511             k = mys[0];
00512             l = mys[1];
00513             if (srfip1_1.idrx != 0) {
00514                 draws_(&m[(i__ + k * m_dim2 << 1) + 1], &m[(i__ + k * m_dim2 
00515                         << 1) + 2], &m[(i__ + l * m_dim2 << 1) + 1], &m[(i__ 
00516                         + l * m_dim2 << 1) + 2], &c__1, &c__1);
00517             }
00518             if (srfblk_1.ndrz != 0 && ii != ni) {
00519 /* Computing MIN */
00520                 i__5 = i__, i__4 = i__ + li;
00521                 i__2 = min(i__5,i__4);
00522                 ctcell_(&z__[z_offset], &mmxx, &nnxx, &nnyy, &m[m_offset], &
00523                         i__2, &k);
00524             }
00525 L320:
00526             i__2 = nypass;
00527             for (jpass = 1; jpass <= i__2; ++jpass) {
00528                 lj = myj[jpass - 1];
00529                 mj = mys[jpass - 1] - lj;
00530                 nj = (i__5 = mj - myf[jpass - 1], abs(i__5));
00531                 i__5 = nj;
00532                 for (jj = 1; jj <= i__5; ++jj) {
00533                     j = mj + jj * lj;
00534                     jplj = j + lj;
00535                     if (srfip1_1.idrx != 0 && jj != nj) {
00536                         draws_(&m[(i__ + j * m_dim2 << 1) + 1], &m[(i__ + j * 
00537                                 m_dim2 << 1) + 2], &m[(i__ + jplj * m_dim2 << 
00538                                 1) + 1], &m[(i__ + jplj * m_dim2 << 1) + 2], &
00539                                 c__1, &c__1);
00540                     }
00541                     if (i__ != mxf[0] && srfip1_1.idry != 0) {
00542                         draws_(&m[(ipli + j * m_dim2 << 1) + 1], &m[(ipli + j 
00543                                 * m_dim2 << 1) + 2], &m[(i__ + j * m_dim2 << 
00544                                 1) + 1], &m[(i__ + j * m_dim2 << 1) + 2], &
00545                                 c__1, &c__1);
00546                     }
00547                     if (srfblk_1.ndrz != 0 && jj != nj && ii != nnxx) {
00548 /* Computing MIN */
00549                         i__6 = i__, i__7 = i__ + li;
00550                         i__4 = min(i__6,i__7);
00551 /* Computing MIN */
00552                         i__9 = j, i__10 = j + lj;
00553                         i__8 = min(i__9,i__10);
00554                         ctcell_(&z__[z_offset], &mmxx, &nnxx, &nnyy, &m[
00555                                 m_offset], &i__4, &i__8);
00556                     }
00557 /* L330: */
00558                 }
00559 /* L340: */
00560             }
00561 /* L350: */
00562         }
00563         goto L430;
00564 L360:
00565         if (srfip1_1.iskirt != 0 || nxpass != 1) {
00566             goto L380;
00567         }
00568         j = mys[0];
00569         i__3 = nnxx;
00570         for (i__ = 2; i__ <= i__3; ++i__) {
00571             draws_(&m[(i__ - 1 + j * m_dim2 << 1) + 1], &m[(i__ - 1 + j * 
00572                     m_dim2 << 1) + 2], &m[(i__ + j * m_dim2 << 1) + 1], &m[(
00573                     i__ + j * m_dim2 << 1) + 2], &c__0, &c__1);
00574 /* L370: */
00575         }
00576 L380:
00577         i__3 = nnyy;
00578         for (jj = 1; jj <= i__3; ++jj) {
00579             j = mj + jj * lj;
00580             jplj = j + lj;
00581             if (nxpass == 1) {
00582                 goto L390;
00583             }
00584             k = mxs[0];
00585             l = mxs[1];
00586             if (srfip1_1.idry != 0) {
00587                 draws_(&m[(k + j * m_dim2 << 1) + 1], &m[(k + j * m_dim2 << 1)
00588                          + 2], &m[(l + j * m_dim2 << 1) + 1], &m[(l + j * 
00589                         m_dim2 << 1) + 2], &c__1, &c__1);
00590             }
00591             if (srfblk_1.ndrz != 0 && jj != nj) {
00592 /* Computing MIN */
00593                 i__5 = j, i__4 = j + lj;
00594                 i__2 = min(i__5,i__4);
00595                 ctcell_(&z__[z_offset], &mmxx, &nnxx, &nnyy, &m[m_offset], &k,
00596                          &i__2);
00597             }
00598 L390:
00599             i__2 = nxpass;
00600             for (ipass = 1; ipass <= i__2; ++ipass) {
00601                 li = mxj[ipass - 1];
00602                 mi = mxs[ipass - 1] - li;
00603                 ni = (i__5 = mi - mxf[ipass - 1], abs(i__5));
00604                 i__5 = ni;
00605                 for (ii = 1; ii <= i__5; ++ii) {
00606                     i__ = mi + ii * li;
00607                     ipli = i__ + li;
00608                     if (srfip1_1.idry != 0 && ii != ni) {
00609                         draws_(&m[(i__ + j * m_dim2 << 1) + 1], &m[(i__ + j * 
00610                                 m_dim2 << 1) + 2], &m[(ipli + j * m_dim2 << 1)
00611                                  + 1], &m[(ipli + j * m_dim2 << 1) + 2], &
00612                                 c__1, &c__1);
00613                     }
00614                     if (j != myf[0] && srfip1_1.idrx != 0) {
00615                         draws_(&m[(i__ + jplj * m_dim2 << 1) + 1], &m[(i__ + 
00616                                 jplj * m_dim2 << 1) + 2], &m[(i__ + j * 
00617                                 m_dim2 << 1) + 1], &m[(i__ + j * m_dim2 << 1) 
00618                                 + 2], &c__1, &c__1);
00619                     }
00620                     if (srfblk_1.ndrz != 0 && ii != ni && jj != nnyy) {
00621 /* Computing MIN */
00622                         i__6 = i__, i__7 = i__ + li;
00623                         i__4 = min(i__6,i__7);
00624 /* Computing MIN */
00625                         i__9 = j, i__10 = j + lj;
00626                         i__8 = min(i__9,i__10);
00627                         ctcell_(&z__[z_offset], &mmxx, &nnxx, &nnyy, &m[
00628                                 m_offset], &i__4, &i__8);
00629                     }
00630 /* L400: */
00631                 }
00632 /* L410: */
00633             }
00634 /* L420: */
00635         }
00636 L430:
00637         if (srfip1_1.iskirt == 0) {
00638             goto L520;
00639         }
00640 
00641 /* FIX UP IF SKIRT IS USED WITH LINES ONE WAY. */
00642 
00643         if (srfip1_1.idrx != 0) {
00644             goto L460;
00645         }
00646         i__3 = nxpass;
00647         for (ipass = 1; ipass <= i__3; ++ipass) {
00648             if (nxpass == 2) {
00649                 if__ = (ipass - 1) * (nnxx - 1) + 1;
00650             }
00651             i__2 = nnyy;
00652             for (j = 2; j <= i__2; ++j) {
00653                 draws_(&m[(if__ + (j - 1) * m_dim2 << 1) + 1], &m[(if__ + (j 
00654                         - 1) * m_dim2 << 1) + 2], &m[(if__ + j * m_dim2 << 1) 
00655                         + 1], &m[(if__ + j * m_dim2 << 1) + 2], &c__1, &c__0);
00656 /* L440: */
00657             }
00658 /* L450: */
00659         }
00660 L460:
00661         if (srfip1_1.idry != 0) {
00662             goto L520;
00663         }
00664         i__3 = nypass;
00665         for (jpass = 1; jpass <= i__3; ++jpass) {
00666             if (nypass == 2) {
00667                 jf = (jpass - 1) * (nnyy - 1) + 1;
00668             }
00669             i__2 = nnxx;
00670             for (i__ = 2; i__ <= i__2; ++i__) {
00671                 draws_(&m[(i__ - 1 + jf * m_dim2 << 1) + 1], &m[(i__ - 1 + jf 
00672                         * m_dim2 << 1) + 2], &m[(i__ + jf * m_dim2 << 1) + 1],
00673                          &m[(i__ + jf * m_dim2 << 1) + 2], &c__1, &c__0);
00674 /* L470: */
00675             }
00676 /* L480: */
00677         }
00678         goto L520;
00679 
00680 /* ALL VISIBLE IF VIEWED FROM DIRECTLY ABOVE OR BELOW. */
00681 
00682 L490:
00683         if (srfblk_1.nupper > 0 && s[3] < s[6]) {
00684             goto L520;
00685         }
00686         if (srfblk_1.nupper < 0 && s[3] > s[6]) {
00687             goto L520;
00688         }
00689         srfblk_1.nupper = 1;
00690         if (s[3] < s[6]) {
00691             srfblk_1.nupper = -1;
00692         }
00693         i__3 = nnxx;
00694         for (i__ = 1; i__ <= i__3; ++i__) {
00695             i__2 = nnyy;
00696             for (j = 1; j <= i__2; ++j) {
00697                 if (srfip1_1.idrx != 0 && j != nnyy) {
00698                     draws_(&m[(i__ + j * m_dim2 << 1) + 1], &m[(i__ + j * 
00699                             m_dim2 << 1) + 2], &m[(i__ + (j + 1) * m_dim2 << 
00700                             1) + 1], &m[(i__ + (j + 1) * m_dim2 << 1) + 2], &
00701                             c__1, &c__0);
00702                 }
00703                 if (srfip1_1.idry != 0 && i__ != nnxx) {
00704                     draws_(&m[(i__ + j * m_dim2 << 1) + 1], &m[(i__ + j * 
00705                             m_dim2 << 1) + 2], &m[(i__ + 1 + j * m_dim2 << 1) 
00706                             + 1], &m[(i__ + 1 + j * m_dim2 << 1) + 2], &c__1, 
00707                             &c__0);
00708                 }
00709                 if (srfip1_1.idrz != 0 && i__ != nnxx && j != nnyy) {
00710                     ctcell_(&z__[z_offset], &mmxx, &nnxx, &nnyy, &m[m_offset],
00711                              &i__, &j);
00712                 }
00713 /* L500: */
00714             }
00715 /* L510: */
00716         }
00717 L520:
00718         if (ster == 0.f) {
00719             goto L560;
00720         }
00721         if (srfip1_1.istp < 0) {
00722             goto L540;
00723         } else if (srfip1_1.istp == 0) {
00724             goto L530;
00725         } else {
00726             goto L550;
00727         }
00728 L530:
00729         frame_();
00730 L540:
00731         frame_();
00732         goto L570;
00733 L550:
00734         if (ipic != 2) {
00735             goto L570;
00736         }
00737 L560:
00738         if (srfip1_1.ifr > 0) {
00739             frame_();
00740         }
00741 L570:
00742         ;
00743     }
00744     return 0;
00745 } /* srface_ */
00746 
00747 
00748 
00749 
00750 /* Subroutine */ int srfpl_(integer *n, real *px, real *py)
00751 {
00752     extern /* Subroutine */ int line_(real *, real *, real *, real *);
00753 
00754     /* Parameter adjustments */
00755     --py;
00756     --px;
00757 
00758     /* Function Body */
00759     line_(&px[1], &py[1], &px[2], &py[2]);
00760     return 0;
00761 } /* srfpl_ */
00762 
00763 
00764 
00765 
00766 /* Subroutine */ int clset_(real *z__, integer *mx, integer *nx, integer *ny, 
00767         real *chi, real *clo, real *cinc, integer *nla, integer *nlm, real *
00768         cl, integer *ncl, integer *icnst, integer *ioffp, real *spval, real *
00769         bigest)
00770 {
00771     /* Initialized data */
00772 
00773     static integer kk = 0;
00774 
00775     /* System generated locals */
00776     integer z_dim1, z_offset, i__1, i__2;
00777     real r__1;
00778 
00779     /* Builtin functions */
00780     double r_lg10(real *), pow_ri(real *, integer *), r_int(real *);
00781 
00782     /* Local variables */
00783     static real fanc, crat;
00784     static integer i__, j, k;
00785     static real p, cc, ha, glo;
00786 
00787 /* cc      DIMENSION       Z(MX,NY)   ,CL(NLM) */
00788     /* Parameter adjustments */
00789     z_dim1 = *mx;
00790     z_offset = z_dim1 + 1;
00791     z__ -= z_offset;
00792     --cl;
00793 
00794     /* Function Body */
00795 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00796  */
00797 /* CLSET PUTS THE VALUS OF THE CONTOUR LEVELS IN CL */
00798 
00799     *icnst = 0;
00800     glo = *clo;
00801     ha = *chi;
00802     fanc = *cinc;
00803     crat = (real) (*nla);
00804     if ((r__1 = ha - glo) < 0.f) {
00805         goto L10;
00806     } else if (r__1 == 0) {
00807         goto L20;
00808     } else {
00809         goto L50;
00810     }
00811 L10:
00812     glo = ha;
00813     ha = *clo;
00814     goto L50;
00815 L20:
00816     glo = *bigest;
00817     ha = -glo;
00818     i__1 = *ny;
00819     for (j = 1; j <= i__1; ++j) {
00820         i__2 = *nx;
00821         for (i__ = 1; i__ <= i__2; ++i__) {
00822             if (*ioffp == 1 && z__[i__ + j * z_dim1] == *spval) {
00823                 goto L30;
00824             }
00825 /* Computing MIN */
00826             r__1 = z__[i__ + j * z_dim1];
00827             glo = dmin(r__1,glo);
00828 /* Computing MAX */
00829             r__1 = z__[i__ + j * z_dim1];
00830             ha = dmax(r__1,ha);
00831 L30:
00832             ;
00833         }
00834 /* L40: */
00835     }
00836 L50:
00837     if (fanc < 0.f) {
00838         goto L60;
00839     } else if (fanc == 0) {
00840         goto L70;
00841     } else {
00842         goto L90;
00843     }
00844 L60:
00845     crat = -fanc;
00846 L70:
00847     fanc = (ha - glo) / crat;
00848     if (fanc <= 0.f) {
00849         goto L140;
00850     } else {
00851         goto L80;
00852     }
00853 L80:
00854     i__1 = (integer) (r_lg10(&fanc) + 500.f) - 500;
00855     p = pow_ri(&c_b128, &i__1);
00856     r__1 = fanc / p;
00857     fanc = r_int(&r__1) * p;
00858 L90:
00859     if (*chi - *clo != 0.f) {
00860         goto L110;
00861     } else {
00862         goto L100;
00863     }
00864 L100:
00865     r__1 = glo / fanc;
00866     glo = r_int(&r__1) * fanc;
00867     r__1 = ha / fanc;
00868     ha = r_int(&r__1) * fanc;
00869 L110:
00870     i__1 = *nlm;
00871     for (k = 1; k <= i__1; ++k) {
00872         cc = glo + (real) (k - 1) * fanc;
00873         if (cc > ha) {
00874             goto L130;
00875         }
00876         kk = k;
00877         cl[k] = cc;
00878 /* L120: */
00879     }
00880 L130:
00881     *ncl = kk;
00882     return 0;
00883 L140:
00884     *icnst = 1;
00885     return 0;
00886 } /* clset_ */
00887 
00888 
00889 
00890 
00891 /* Subroutine */ int ctcell_(real *z__, integer *mx, integer *nx, integer *ny,
00892          integer *m, integer *i0, integer *j0)
00893 {
00894     /* Initialized data */
00895 
00896     static integer idub = 0;
00897 
00898     /* System generated locals */
00899     integer z_dim1, z_offset, m_dim2, m_offset, i__1, i__2;
00900     real r__1;
00901 
00902     /* Builtin functions */
00903     double r_sign(real *, real *);
00904 
00905     /* Local variables */
00906     static integer jump, k;
00907     extern /* Subroutine */ int color_(integer *), draws_(integer *, integer *
00908             , integer *, integer *, integer *, integer *);
00909     static integer i1, j1;
00910     static real h1, h2, h3, h4;
00911     static integer k1, k2, k3, k4;
00912     static real ra, rb, cv;
00913     static logical lcolor;
00914     static integer i1p1, j1p1, mua, mva, mub, mvb;
00915 
00916 
00917 /* CC      DIMENSION       Z(MX,NY)   ,M(2,NX,NY) */
00918 
00919     /* Parameter adjustments */
00920     m_dim2 = *nx;
00921     m_offset = (m_dim2 + 1 << 1) + 1;
00922     m -= m_offset;
00923     z_dim1 = *mx;
00924     z_offset = z_dim1 + 1;
00925     z__ -= z_offset;
00926 
00927     /* Function Body */
00928 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00929  */
00930 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00931  */
00932     i1 = *i0;
00933     i1p1 = i1 + 1;
00934     j1 = *j0;
00935     j1p1 = j1 + 1;
00936     h1 = z__[i1 + j1 * z_dim1];
00937     h2 = z__[i1 + j1p1 * z_dim1];
00938     h3 = z__[i1p1 + j1p1 * z_dim1];
00939     h4 = z__[i1p1 + j1 * z_dim1];
00940     if (srfblk_1.ioffp != 1) {
00941         goto L10;
00942     }
00943     if (h1 == srfblk_1.spval || h2 == srfblk_1.spval || h3 == srfblk_1.spval 
00944             || h4 == srfblk_1.spval) {
00945         return 0;
00946     }
00947 L10:
00948 /* Computing MIN */
00949     r__1 = min(h1,h2), r__1 = min(r__1,h3);
00950     if (dmin(r__1,h4) > srfblk_1.cl[srfblk_1.ncl - 1]) {
00951         return 0;
00952     }
00953 
00954     lcolor = FALSE_;
00955     i__1 = srfblk_1.ncl;
00956     for (k = 1; k <= i__1; ++k) {
00957 
00958 /* FOR EACH CONTOUR LEVEL, DESIDE WHICH OF THE 16 BASIC SIT- */
00959 /* UATIONS EXISTS, THEN INTERPOLATE IN TWO-SPACE TO FIND THE */
00960 /* END POINTS OF THE CONTOUR LINE SEGMENT WITHIN THIS CELL. */
00961 
00962         cv = srfblk_1.cl[k - 1];
00963         r__1 = h1 - cv;
00964         k1 = ((integer) r_sign(&c_b3, &r__1) + 1) / 2;
00965         r__1 = h2 - cv;
00966         k2 = ((integer) r_sign(&c_b3, &r__1) + 1) / 2;
00967         r__1 = h3 - cv;
00968         k3 = ((integer) r_sign(&c_b3, &r__1) + 1) / 2;
00969         r__1 = h4 - cv;
00970         k4 = ((integer) r_sign(&c_b3, &r__1) + 1) / 2;
00971         jump = k1 + 1 + (k2 << 1) + (k3 << 2) + (k4 << 3);
00972 
00973 /*  17/Apr/91:  plot contours in different colors */
00974 
00975         if (jump > 1 && jump < 16) {
00976             i__2 = k % 6 + 2;
00977             color_(&i__2);
00978         }
00979 
00980         switch (jump) {
00981             case 1:  goto L120;
00982             case 2:  goto L30;
00983             case 3:  goto L50;
00984             case 4:  goto L60;
00985             case 5:  goto L70;
00986             case 6:  goto L20;
00987             case 7:  goto L80;
00988             case 8:  goto L90;
00989             case 9:  goto L90;
00990             case 10:  goto L80;
00991             case 11:  goto L40;
00992             case 12:  goto L70;
00993             case 13:  goto L60;
00994             case 14:  goto L50;
00995             case 15:  goto L30;
00996             case 16:  goto L110;
00997         }
00998 
00999 L20:
01000         idub = 1;
01001 L30:
01002         ra = (h1 - cv) / (h1 - h2);
01003         mua = (real) m[(i1 + j1 * m_dim2 << 1) + 1] + ra * (real) (m[(i1 + 
01004                 j1p1 * m_dim2 << 1) + 1] - m[(i1 + j1 * m_dim2 << 1) + 1]);
01005         mva = (real) m[(i1 + j1 * m_dim2 << 1) + 2] + ra * (real) (m[(i1 + 
01006                 j1p1 * m_dim2 << 1) + 2] - m[(i1 + j1 * m_dim2 << 1) + 2]);
01007         rb = (h1 - cv) / (h1 - h4);
01008         mub = (real) m[(i1 + j1 * m_dim2 << 1) + 1] + rb * (real) (m[(i1p1 + 
01009                 j1 * m_dim2 << 1) + 1] - m[(i1 + j1 * m_dim2 << 1) + 1]);
01010         mvb = (real) m[(i1 + j1 * m_dim2 << 1) + 2] + rb * (real) (m[(i1p1 + 
01011                 j1 * m_dim2 << 1) + 2] - m[(i1 + j1 * m_dim2 << 1) + 2]);
01012         goto L100;
01013 L40:
01014         idub = -1;
01015 L50:
01016         ra = (h2 - cv) / (h2 - h1);
01017         mua = (real) m[(i1 + j1p1 * m_dim2 << 1) + 1] + ra * (real) (m[(i1 + 
01018                 j1 * m_dim2 << 1) + 1] - m[(i1 + j1p1 * m_dim2 << 1) + 1]);
01019         mva = (real) m[(i1 + j1p1 * m_dim2 << 1) + 2] + ra * (real) (m[(i1 + 
01020                 j1 * m_dim2 << 1) + 2] - m[(i1 + j1p1 * m_dim2 << 1) + 2]);
01021         rb = (h2 - cv) / (h2 - h3);
01022         mub = (real) m[(i1 + j1p1 * m_dim2 << 1) + 1] + rb * (real) (m[(i1p1 
01023                 + j1p1 * m_dim2 << 1) + 1] - m[(i1 + j1p1 * m_dim2 << 1) + 1])
01024                 ;
01025         mvb = (real) m[(i1 + j1p1 * m_dim2 << 1) + 2] + rb * (real) (m[(i1p1 
01026                 + j1p1 * m_dim2 << 1) + 2] - m[(i1 + j1p1 * m_dim2 << 1) + 2])
01027                 ;
01028         goto L100;
01029 L60:
01030         ra = (h2 - cv) / (h2 - h3);
01031         mua = (real) m[(i1 + j1p1 * m_dim2 << 1) + 1] + ra * (real) (m[(i1p1 
01032                 + j1p1 * m_dim2 << 1) + 1] - m[(i1 + j1p1 * m_dim2 << 1) + 1])
01033                 ;
01034         mva = (real) m[(i1 + j1p1 * m_dim2 << 1) + 2] + ra * (real) (m[(i1p1 
01035                 + j1p1 * m_dim2 << 1) + 2] - m[(i1 + j1p1 * m_dim2 << 1) + 2])
01036                 ;
01037         rb = (h1 - cv) / (h1 - h4);
01038         mub = (real) m[(i1 + j1 * m_dim2 << 1) + 1] + rb * (real) (m[(i1p1 + 
01039                 j1 * m_dim2 << 1) + 1] - m[(i1 + j1 * m_dim2 << 1) + 1]);
01040         mvb = (real) m[(i1 + j1 * m_dim2 << 1) + 2] + rb * (real) (m[(i1p1 + 
01041                 j1 * m_dim2 << 1) + 2] - m[(i1 + j1 * m_dim2 << 1) + 2]);
01042         goto L100;
01043 L70:
01044         ra = (h3 - cv) / (h3 - h2);
01045         mua = (real) m[(i1p1 + j1p1 * m_dim2 << 1) + 1] + ra * (real) (m[(i1 
01046                 + j1p1 * m_dim2 << 1) + 1] - m[(i1p1 + j1p1 * m_dim2 << 1) + 
01047                 1]);
01048         mva = (real) m[(i1p1 + j1p1 * m_dim2 << 1) + 2] + ra * (real) (m[(i1 
01049                 + j1p1 * m_dim2 << 1) + 2] - m[(i1p1 + j1p1 * m_dim2 << 1) + 
01050                 2]);
01051         rb = (h3 - cv) / (h3 - h4);
01052         mub = (real) m[(i1p1 + j1p1 * m_dim2 << 1) + 1] + rb * (real) (m[(
01053                 i1p1 + j1 * m_dim2 << 1) + 1] - m[(i1p1 + j1p1 * m_dim2 << 1) 
01054                 + 1]);
01055         mvb = (real) m[(i1p1 + j1p1 * m_dim2 << 1) + 2] + rb * (real) (m[(
01056                 i1p1 + j1 * m_dim2 << 1) + 2] - m[(i1p1 + j1p1 * m_dim2 << 1) 
01057                 + 2]);
01058         idub = 0;
01059         goto L100;
01060 L80:
01061         ra = (h2 - cv) / (h2 - h1);
01062         mua = (real) m[(i1 + j1p1 * m_dim2 << 1) + 1] + ra * (real) (m[(i1 + 
01063                 j1 * m_dim2 << 1) + 1] - m[(i1 + j1p1 * m_dim2 << 1) + 1]);
01064         mva = (real) m[(i1 + j1p1 * m_dim2 << 1) + 2] + ra * (real) (m[(i1 + 
01065                 j1 * m_dim2 << 1) + 2] - m[(i1 + j1p1 * m_dim2 << 1) + 2]);
01066         rb = (h3 - cv) / (h3 - h4);
01067         mub = (real) m[(i1p1 + j1p1 * m_dim2 << 1) + 1] + rb * (real) (m[(
01068                 i1p1 + j1 * m_dim2 << 1) + 1] - m[(i1p1 + j1p1 * m_dim2 << 1) 
01069                 + 1]);
01070         mvb = (real) m[(i1p1 + j1p1 * m_dim2 << 1) + 2] + rb * (real) (m[(
01071                 i1p1 + j1 * m_dim2 << 1) + 2] - m[(i1p1 + j1p1 * m_dim2 << 1) 
01072                 + 2]);
01073         goto L100;
01074 L90:
01075         ra = (h4 - cv) / (h4 - h1);
01076         mua = (real) m[(i1p1 + j1 * m_dim2 << 1) + 1] + ra * (real) (m[(i1 + 
01077                 j1 * m_dim2 << 1) + 1] - m[(i1p1 + j1 * m_dim2 << 1) + 1]);
01078         mva = (real) m[(i1p1 + j1 * m_dim2 << 1) + 2] + ra * (real) (m[(i1 + 
01079                 j1 * m_dim2 << 1) + 2] - m[(i1p1 + j1 * m_dim2 << 1) + 2]);
01080         rb = (h4 - cv) / (h4 - h3);
01081         mub = (real) m[(i1p1 + j1 * m_dim2 << 1) + 1] + rb * (real) (m[(i1p1 
01082                 + j1p1 * m_dim2 << 1) + 1] - m[(i1p1 + j1 * m_dim2 << 1) + 1])
01083                 ;
01084         mvb = (real) m[(i1p1 + j1 * m_dim2 << 1) + 2] + rb * (real) (m[(i1p1 
01085                 + j1p1 * m_dim2 << 1) + 2] - m[(i1p1 + j1 * m_dim2 << 1) + 2])
01086                 ;
01087         idub = 0;
01088 L100:
01089         draws_(&mua, &mva, &mub, &mvb, &c__1, &c__0);
01090         lcolor = TRUE_;
01091         if (idub < 0) {
01092             goto L90;
01093         } else if (idub == 0) {
01094             goto L110;
01095         } else {
01096             goto L70;
01097         }
01098 L110:
01099         ;
01100     }
01101 
01102 L120:
01103     if (lcolor) {
01104         color_(&c__1);
01105     }
01106     return 0;
01107 } /* ctcell_ */
01108 
01109 
01110 
01111 
01112 /* Subroutine */ int draws_(integer *mx1, integer *my1, integer *mx2, integer 
01113         *my2, integer *idraw, integer *imark)
01114 {
01115     /* Initialized data */
01116 
01117     static real steep = 5.f;
01118     static integer mx = 0;
01119     static integer my = 0;
01120 
01121     /* System generated locals */
01122     integer i__1, i__2;
01123 
01124     /* Local variables */
01125     static integer nx1p1, k, ltemp;
01126     extern /* Subroutine */ int srfpl_(integer *, real *, real *);
01127     static real dy;
01128     static integer nx1, ny1, nx2, ny2;
01129     static real pxs[2], pys[2], fny1;
01130     static logical vis1, vis2;
01131     static integer mmx1, mmy1, mmx2, mmy2;
01132 
01133 
01134 /* THIS ROUTINE DRAWS THE VISIBLE PART OF THE LINE CONNECTING */
01135 /* (MX1,MY1) AND (MX2,MY2).  IF IDRAW .NE. 0, THE LINE IS DRAWN. */
01136 /* IF IMARK .NE. 0, THE VISIBILITY ARRAY IS MARKED. */
01137 
01138 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01139  */
01140 /* MAKE LINE LEFT TO RIGHT. */
01141 
01142     mmx1 = *mx1;
01143     mmy1 = *my1;
01144     mmx2 = *mx2;
01145     mmy2 = *my2;
01146     if (mmx1 == srfblk_1.nspval || mmx2 == srfblk_1.nspval) {
01147         return 0;
01148     }
01149     if (mmx1 > mmx2) {
01150         goto L10;
01151     }
01152     nx1 = mmx1;
01153     ny1 = mmy1;
01154     nx2 = mmx2;
01155     ny2 = mmy2;
01156     goto L20;
01157 L10:
01158     nx1 = mmx2;
01159     ny1 = mmy2;
01160     nx2 = mmx1;
01161     ny2 = mmy1;
01162 L20:
01163     if (srfblk_1.nupper < 0) {
01164         goto L180;
01165     }
01166 
01167 /* CHECK UPPER VISIBILITY. */
01168 
01169     vis1 = ny1 >= srfblk_1.limu[nx1 - 1] - 1;
01170     vis2 = ny2 >= srfblk_1.limu[nx2 - 1] - 1;
01171 
01172 /* VIS1 AND VIS2 TRUE MEANS VISIBLE. */
01173 
01174     if (vis1 && vis2) {
01175         goto L120;
01176     }
01177 
01178 /* VIS1 AND VIS2 FALSE MEANS INVISIBLE. */
01179 
01180     if (! (vis1 || vis2)) {
01181         goto L180;
01182     }
01183 
01184 /* FIND CHANGE POINT. */
01185 
01186     if (nx1 == nx2) {
01187         goto L110;
01188     }
01189     dy = (real) (ny2 - ny1) / (real) (nx2 - nx1);
01190     nx1p1 = nx1 + 1;
01191     fny1 = (real) ny1;
01192     if (vis1) {
01193         goto L60;
01194     }
01195     i__1 = nx2;
01196     for (k = nx1p1; k <= i__1; ++k) {
01197         mx = k;
01198         my = fny1 + (real) (k - nx1) * dy;
01199         if (my > srfblk_1.limu[k - 1]) {
01200             goto L40;
01201         }
01202 /* L30: */
01203     }
01204 L40:
01205     if (dabs(dy) >= steep) {
01206         goto L90;
01207     }
01208 L50:
01209     nx1 = mx;
01210     ny1 = my;
01211     goto L120;
01212 L60:
01213     i__1 = nx2;
01214     for (k = nx1p1; k <= i__1; ++k) {
01215         mx = k;
01216         my = fny1 + (real) (k - nx1) * dy;
01217         if (my <= srfblk_1.limu[k - 1]) {
01218             goto L80;
01219         }
01220 /* L70: */
01221     }
01222 L80:
01223     if (dabs(dy) >= steep) {
01224         goto L100;
01225     }
01226     nx2 = mx - 1;
01227     ny2 = my;
01228     goto L120;
01229 L90:
01230     if (srfblk_1.limu[mx - 1] == 0) {
01231         goto L50;
01232     }
01233     nx1 = mx;
01234     ny1 = srfblk_1.limu[nx1 - 1];
01235     goto L120;
01236 L100:
01237     nx2 = mx - 1;
01238     ny2 = srfblk_1.limu[nx2 - 1];
01239     goto L120;
01240 L110:
01241     if (vis1) {
01242 /* Computing MIN */
01243         i__1 = srfblk_1.limu[nx1 - 1], i__2 = srfblk_1.limu[nx2 - 1];
01244         ny2 = min(i__1,i__2);
01245     }
01246     if (vis2) {
01247 /* Computing MIN */
01248         i__1 = srfblk_1.limu[nx1 - 1], i__2 = srfblk_1.limu[nx2 - 1];
01249         ny1 = min(i__1,i__2);
01250     }
01251 L120:
01252     if (*idraw == 0) {
01253         goto L150;
01254     }
01255 
01256 /* DRAW VISIBLE PART OF LINE. */
01257 
01258     if (srfblk_1.irot != 0) {
01259         goto L130;
01260     } else {
01261         goto L140;
01262     }
01263 L130:
01264     pxs[0] = (real) ny1;
01265     pxs[1] = (real) ny2;
01266     pys[0] = (real) (1024 - nx1);
01267     pys[1] = (real) (1024 - nx2);
01268     srfpl_(&c__2, pxs, pys);
01269     goto L150;
01270 L140:
01271     pxs[0] = (real) nx1;
01272     pxs[1] = (real) nx2;
01273     pys[0] = (real) ny1;
01274     pys[1] = (real) ny2;
01275     srfpl_(&c__2, pxs, pys);
01276 L150:
01277     if (*imark == 0) {
01278         goto L180;
01279     }
01280     if (nx1 == nx2) {
01281         goto L170;
01282     }
01283     dy = (real) (ny2 - ny1) / (real) (nx2 - nx1);
01284     fny1 = (real) ny1;
01285     i__1 = nx2;
01286     for (k = nx1; k <= i__1; ++k) {
01287         ltemp = fny1 + (real) (k - nx1) * dy;
01288         if (ltemp > srfblk_1.limu[k - 1]) {
01289             srfblk_1.limu[k - 1] = ltemp;
01290         }
01291 /* L160: */
01292     }
01293     goto L180;
01294 L170:
01295     ltemp = max(ny1,ny2);
01296     if (ltemp > srfblk_1.limu[nx1 - 1]) {
01297         srfblk_1.limu[nx1 - 1] = ltemp;
01298     }
01299 L180:
01300     if (srfblk_1.nupper <= 0) {
01301         goto L190;
01302     } else {
01303         goto L370;
01304     }
01305 
01306 /* SAME IDEA AS ABOVE, BUT FOR LOWER SIDE. */
01307 
01308 L190:
01309     if (mmx1 > mmx2) {
01310         goto L200;
01311     }
01312     nx1 = mmx1;
01313     ny1 = mmy1;
01314     nx2 = mmx2;
01315     ny2 = mmy2;
01316     goto L210;
01317 L200:
01318     nx1 = mmx2;
01319     ny1 = mmy2;
01320     nx2 = mmx1;
01321     ny2 = mmy1;
01322 L210:
01323     vis1 = ny1 <= srfblk_1.liml[nx1 - 1] + 1;
01324     vis2 = ny2 <= srfblk_1.liml[nx2 - 1] + 1;
01325     if (vis1 && vis2) {
01326         goto L310;
01327     }
01328     if (! (vis1 || vis2)) {
01329         goto L370;
01330     }
01331     if (nx1 == nx2) {
01332         goto L300;
01333     }
01334     dy = (real) (ny2 - ny1) / (real) (nx2 - nx1);
01335     nx1p1 = nx1 + 1;
01336     fny1 = (real) ny1;
01337     if (vis1) {
01338         goto L250;
01339     }
01340     i__1 = nx2;
01341     for (k = nx1p1; k <= i__1; ++k) {
01342         mx = k;
01343         my = fny1 + (real) (k - nx1) * dy;
01344         if (my < srfblk_1.liml[k - 1]) {
01345             goto L230;
01346         }
01347 /* L220: */
01348     }
01349 L230:
01350     if (dabs(dy) >= steep) {
01351         goto L280;
01352     }
01353 L240:
01354     nx1 = mx;
01355     ny1 = my;
01356     goto L310;
01357 L250:
01358     i__1 = nx2;
01359     for (k = nx1p1; k <= i__1; ++k) {
01360         mx = k;
01361         my = fny1 + (real) (k - nx1) * dy;
01362         if (my >= srfblk_1.liml[k - 1]) {
01363             goto L270;
01364         }
01365 /* L260: */
01366     }
01367 L270:
01368     if (dabs(dy) >= steep) {
01369         goto L290;
01370     }
01371     nx2 = mx - 1;
01372     ny2 = my;
01373     goto L310;
01374 L280:
01375     if (srfblk_1.liml[mx - 1] == 1024) {
01376         goto L240;
01377     }
01378     nx1 = mx;
01379     ny1 = srfblk_1.liml[nx1 - 1];
01380     goto L310;
01381 L290:
01382     nx2 = mx - 1;
01383     ny2 = srfblk_1.liml[nx2 - 1];
01384     goto L310;
01385 L300:
01386     if (vis1) {
01387 /* Computing MAX */
01388         i__1 = srfblk_1.liml[nx1 - 1], i__2 = srfblk_1.liml[nx2 - 1];
01389         ny2 = max(i__1,i__2);
01390     }
01391     if (vis2) {
01392 /* Computing MAX */
01393         i__1 = srfblk_1.liml[nx1 - 1], i__2 = srfblk_1.liml[nx2 - 1];
01394         ny1 = max(i__1,i__2);
01395     }
01396 L310:
01397     if (*idraw == 0) {
01398         goto L340;
01399     }
01400     if (srfblk_1.irot != 0) {
01401         goto L320;
01402     } else {
01403         goto L330;
01404     }
01405 L320:
01406     pxs[0] = (real) ny1;
01407     pxs[1] = (real) ny2;
01408     pys[0] = (real) (1024 - nx1);
01409     pys[1] = (real) (1024 - nx2);
01410     srfpl_(&c__2, pxs, pys);
01411     goto L340;
01412 L330:
01413     pxs[0] = (real) nx1;
01414     pxs[1] = (real) nx2;
01415     pys[0] = (real) ny1;
01416     pys[1] = (real) ny2;
01417     srfpl_(&c__2, pxs, pys);
01418 L340:
01419     if (*imark == 0) {
01420         goto L370;
01421     }
01422     if (nx1 == nx2) {
01423         goto L360;
01424     }
01425     dy = (real) (ny2 - ny1) / (real) (nx2 - nx1);
01426     fny1 = (real) ny1;
01427     i__1 = nx2;
01428     for (k = nx1; k <= i__1; ++k) {
01429         ltemp = fny1 + (real) (k - nx1) * dy;
01430         if (ltemp < srfblk_1.liml[k - 1]) {
01431             srfblk_1.liml[k - 1] = ltemp;
01432         }
01433 /* L350: */
01434     }
01435     return 0;
01436 L360:
01437     ltemp = min(ny1,ny2);
01438     if (ltemp < srfblk_1.liml[nx1 - 1]) {
01439         srfblk_1.liml[nx1 - 1] = ltemp;
01440     }
01441 L370:
01442     return 0;
01443 } /* draws_ */
01444 
01445 
01446 
01447 
01448 /* Subroutine */ int setr_(real *xmin, real *xmax, real *ymin, real *ymax, 
01449         real *zmin, real *zmax, real *r0)
01450 {
01451     /* System generated locals */
01452     real r__1, r__2, r__3;
01453 
01454     /* Local variables */
01455     static real yeye, xeye, zeye, alpha;
01456     extern /* Subroutine */ int trn32s_(real *, real *, real *, real *, real *
01457             , real *, integer *);
01458     static real dummy, dummie, xat, yat, zat, umn, vmn, xmn, ymn, zmn, umx, 
01459             vmx, xmx, ymx, zmx;
01460 
01461 
01462 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01463  */
01464     if (*r0 <= 0.f) {
01465         goto L10;
01466     } else {
01467         goto L20;
01468     }
01469 L10:
01470     srfblk_1.nrswt = 0;
01471     return 0;
01472 L20:
01473     srfblk_1.nrswt = 1;
01474     pwrz1s_1.xxmin = *xmin;
01475     pwrz1s_1.xxmax = *xmax;
01476     pwrz1s_1.yymin = *ymin;
01477     pwrz1s_1.yymax = *ymax;
01478     pwrz1s_1.zzmin = *zmin;
01479     pwrz1s_1.zzmax = *zmax;
01480     srfblk_1.rzero = *r0;
01481     srfblk_1.ll = 0;
01482     xat = (pwrz1s_1.xxmax + pwrz1s_1.xxmin) * .5f;
01483     yat = (pwrz1s_1.yymax + pwrz1s_1.yymin) * .5f;
01484     zat = (pwrz1s_1.zzmax + pwrz1s_1.zzmin) * .5f;
01485     alpha = -(pwrz1s_1.yymin - yat) / (pwrz1s_1.xxmin - xat);
01486     yeye = -srfblk_1.rzero / sqrt(alpha * alpha + 1.f);
01487     xeye = yeye * alpha;
01488     yeye += yat;
01489     xeye += xat;
01490     zeye = zat;
01491     trn32s_(&xat, &yat, &zat, &xeye, &yeye, &zeye, &c__0);
01492     xmn = pwrz1s_1.xxmin;
01493     xmx = pwrz1s_1.xxmax;
01494     ymn = pwrz1s_1.yymin;
01495     ymx = pwrz1s_1.yymax;
01496     zmn = pwrz1s_1.zzmin;
01497     zmx = pwrz1s_1.zzmax;
01498     trn32s_(&xmn, &ymn, &zat, &umn, &dummy, &dummie, &c__1);
01499     trn32s_(&xmx, &ymn, &zmn, &dummy, &vmn, &dummie, &c__1);
01500     trn32s_(&xmx, &ymx, &zat, &umx, &dummy, &dummie, &c__1);
01501     trn32s_(&xmx, &ymn, &zmx, &dummy, &vmx, &dummie, &c__1);
01502     srfblk_1.umin = umn;
01503     srfblk_1.umax = umx;
01504     srfblk_1.vmin = vmn;
01505     srfblk_1.vmax = vmx;
01506 /* Computing 2nd power */
01507     r__1 = pwrz1s_1.xxmax - pwrz1s_1.xxmin;
01508 /* Computing 2nd power */
01509     r__2 = pwrz1s_1.yymax - pwrz1s_1.yymin;
01510 /* Computing 2nd power */
01511     r__3 = pwrz1s_1.zzmax - pwrz1s_1.zzmin;
01512     srfblk_1.bigd = sqrt(r__1 * r__1 + r__2 * r__2 + r__3 * r__3) * .5f;
01513     return 0;
01514 } /* setr_ */
01515 
01516 
01517 
01518 
01519 /* Subroutine */ int trn32s_(real *x, real *y, real *z__, real *xt, real *yt, 
01520         real *zt, integer *iflag)
01521 {
01522     /* Initialized data */
01523 
01524     static integer nlu[7] = { 10,10,100,10,10,10,512 };
01525     static integer nru[7] = { 1014,924,1014,1014,1014,512,1014 };
01526     static integer nbv[7] = { 10,50,50,10,10,256,256 };
01527     static integer ntv[7] = { 1014,964,964,1014,1014,758,758 };
01528 
01529     /* Format strings */
01530     static char fmt_60[] = "";
01531     static char fmt_50[] = "";
01532     static char fmt_120[] = "";
01533     static char fmt_100[] = "";
01534     static char fmt_70[] = "";
01535     static char fmt_80[] = "";
01536 
01537     /* System generated locals */
01538     real r__1, r__2, r__3, r__4;
01539 
01540     /* Local variables */
01541     static integer jump, jump2, jump3;
01542     static real d__, q, r__, cosbe, cosga, sinbe, cosal, singa, u0, v0, u1, 
01543             v1, u2, v2, u3, v3, u4, v4, ax, ay, az, dx, ex, ey, ez, dy, dz, 
01544             xx, yy, zz;
01545 
01546     /* Assigned format variables */
01547     static char *jump3_fmt, *jump2_fmt, *jump_fmt;
01548 
01549 
01550 
01551 /* PICTURE CORNER COORDINATES FOR LL=1 */
01552 
01553 
01554 /* PICTURE CORNER COORDINATES FOR LL=2 */
01555 
01556 
01557 /* PICTURE CORNER COORDINATES FOR LL=3 */
01558 
01559 
01560 /* PICTURE CORNER COORDINATES FOR LL=4 */
01561 
01562 
01563 /* PICTURE CORNER COORDINATES FOR LL=5 */
01564 
01565 
01566 /* PICTURE CORNER COORDINATES FOR LL=6 */
01567 
01568 
01569 /* PICTURE CORNER COORDINATES FOR LL=7 */
01570 
01571 /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
01572  */
01573 /* STORE THE PARAMETERS OF THE SET32 CALL FOR USE WHEN */
01574 /* TRN32 IS CALLED. */
01575 
01576     if (*iflag != 0) {
01577         goto L40;
01578     } else {
01579         goto L10;
01580     }
01581 L10:
01582     jump3 = 0;
01583     jump3_fmt = fmt_60;
01584     if (srfblk_1.ioffp == 1) {
01585         jump3 = 1;
01586         jump3_fmt = fmt_50;
01587     }
01588     ax = *x;
01589     ay = *y;
01590     az = *z__;
01591     ex = *xt;
01592     ey = *yt;
01593     ez = *zt;
01594 
01595 /* AS MUCH COMPUTATION AS POSSIBLE IS DONE DURING EXECUTION */
01596 /* THIS ROUTINE WHEN IFLAG=0 BECAUSE CALLS IN THAT MODE ARE INFREQUENT. */
01597 
01598     dx = ax - ex;
01599     dy = ay - ey;
01600     dz = az - ez;
01601     d__ = sqrt(dx * dx + dy * dy + dz * dz);
01602     cosal = dx / d__;
01603     cosbe = dy / d__;
01604     cosga = dz / d__;
01605     singa = sqrt(1.f - cosga * cosga);
01606     jump2 = 0;
01607     jump2_fmt = fmt_120;
01608     if (srfblk_1.ll == 0) {
01609         goto L20;
01610     }
01611     jump2 = 1;
01612     jump2_fmt = fmt_100;
01613     pwrz1s_1.delcrt = (real) (nru[srfblk_1.ll - 1] - nlu[srfblk_1.ll - 1]);
01614     u0 = srfblk_1.umin;
01615     v0 = srfblk_1.vmin;
01616     u1 = (real) nlu[srfblk_1.ll - 1];
01617     v1 = (real) nbv[srfblk_1.ll - 1];
01618     u2 = (real) (nru[srfblk_1.ll - 1] - nlu[srfblk_1.ll - 1]);
01619     v2 = (real) (ntv[srfblk_1.ll - 1] - nbv[srfblk_1.ll - 1]);
01620     u3 = u2 / (srfblk_1.umax - srfblk_1.umin);
01621     v3 = v2 / (srfblk_1.vmax - srfblk_1.vmin);
01622     u4 = (real) nru[srfblk_1.ll - 1];
01623     v4 = (real) ntv[srfblk_1.ll - 1];
01624     if (srfblk_1.nrswt == 0) {
01625         goto L20;
01626     }
01627     u0 = -srfblk_1.bigd;
01628     v0 = -srfblk_1.bigd;
01629     u3 = u2 / (srfblk_1.bigd * 2.f);
01630     v3 = v2 / (srfblk_1.bigd * 2.f);
01631 
01632 /* THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF */
01633 /* THE 2-SPACE.  THE 3-SPACE Z AXIS IS TRANSFORMED INTO THE */
01634 /* 2-SPACE Y AXIS.  IF THE LINE OF SIGHT IS CLOSE TO PARALLEL */
01635 /* TO THE 3-SPACE Z AXIS, THE 3-SPACE Y AXIS IS CHOSEN (IN- */
01636 /* STEAD OF THE 3-SPACE Z AXIS) TO BE TRANSFORMED INTO THE */
01637 /* 2-SPACE Y AXIS. */
01638 
01639 L20:
01640     if (singa < 1e-4f) {
01641         goto L30;
01642     }
01643     r__ = 1.f / singa;
01644     jump = 0;
01645     jump_fmt = fmt_70;
01646     return 0;
01647 L30:
01648     sinbe = sqrt(1.f - cosbe * cosbe);
01649     r__ = 1.f / sinbe;
01650     jump = 1;
01651     jump_fmt = fmt_80;
01652     return 0;
01653 L40:
01654     xx = *x;
01655     yy = *y;
01656     zz = *z__;
01657     switch (jump3) {
01658         case 0: goto L60;
01659         case 1: goto L50;
01660     }
01661 L50:
01662     if (zz == srfblk_1.spval) {
01663         goto L110;
01664     }
01665 L60:
01666     q = d__ / ((xx - ex) * cosal + (yy - ey) * cosbe + (zz - ez) * cosga);
01667     switch (jump) {
01668         case 0: goto L70;
01669         case 1: goto L80;
01670     }
01671 L70:
01672     xx = ((ex + q * (xx - ex) - ax) * cosbe - (ey + q * (yy - ey) - ay) * 
01673             cosal) * r__;
01674     yy = (ez + q * (zz - ez) - az) * r__;
01675     goto L90;
01676 L80:
01677     xx = ((ez + q * (zz - ez) - az) * cosal - (ex + q * (xx - ex) - ax) * 
01678             cosga) * r__;
01679     yy = (ey + q * (yy - ey) - ay) * r__;
01680 L90:
01681     switch (jump2) {
01682         case 0: goto L120;
01683         case 1: goto L100;
01684     }
01685 L100:
01686 /* Computing MIN */
01687 /* Computing MAX */
01688     r__3 = u1, r__4 = u1 + u3 * (srfblk_1.fact * xx - u0);
01689     r__1 = u4, r__2 = dmax(r__3,r__4);
01690     xx = dmin(r__1,r__2);
01691 /* Computing MIN */
01692 /* Computing MAX */
01693     r__3 = v1, r__4 = v1 + v3 * (srfblk_1.fact * yy - v0);
01694     r__1 = v4, r__2 = dmax(r__3,r__4);
01695     yy = dmin(r__1,r__2);
01696     goto L120;
01697 L110:
01698     xx = (real) srfblk_1.nspval;
01699     yy = (real) srfblk_1.nspval;
01700 
01701 L120:
01702     *xt = xx;
01703     *yt = yy;
01704     return 0;
01705 } /* trn32s_ */
01706 
01707 
01708 /* cc      BLOCKDATA SRFABD */
01709 /* Subroutine */ int srfabd_(void)
01710 {
01711 
01712 /*  INITIALIZATION OF INTERNAL PARAMETERS */
01713 
01714 
01715     return 0;
01716 } /* srfabd_ */
01717