import java.io.*;
import java.applet.*;
import java.awt.*;
import java.awt.event.*;
import java.awt.font.*;
import java.net.URL;
import java.text.*;
import java.util.*;
import java.util.Map;
import java.util.Hashtable;
import java.util.jar.Attributes;


public class Isolines {
    private static int NMAX = 5000;  //was (NMAX=10000)

    private static int xscale = 6;   // number of real pixels per "pixel"
    private static int yscale = 6;   // number of real pixels per "pixel"


    // Draws isolines of the given 2D function at the given function values.
    // Each of the contour lines is made up of a number of segments,
    // each of which goes from one edge of a "pixel" to another.
    // The pixels are defined as the edges of the grid defined by {\tt xvec} 
    //  and {\tt yvec}.  (NOTE: a "pixel" is really multiple pixels - scaled)

    // REAL z(nx*ny) ! arranged: (iy-1)*nx + ix
    public void plot_contours_rect(Graphics g, double z[], double xvec[],
                                          double yvec[], int nx, int ny,
                                          double xmin, double xmax, 
                                          double ymin, double ymax,
                                          int ncontours, double contours[]) {
        int i, j, jj;
        int L;
        double xdenom, ydenom;
        double xmin1, xmax1, ymin1, ymax1;
        //double xscale, yscale;
        int jstart, jend;

        double zvalue;
        int n;
        double x[];
        double y[];

        //xscale = p_plotsize_inches(1);
        //yscale = p_plotsize_inches(2);

        // scale the input vectors
        xdenom = xmax - xmin;
        ydenom = ymax - ymin;

        if (xdenom == 0 || ydenom == 0) {
            // Need to auto-scale the vectors here
            xmax1 = xvec[0];
            xmin1 = xvec[0];
            ymax1 = yvec[0];
            ymin1 = yvec[0];
            for (i = 0; i < nx; i++) {
                if (xvec[i] > xmax1) xmax1 = xvec[i];
                if (xvec[i] < xmin1) xmin1 = xvec[i];
            }
            for (i = 0; i < ny; i++) {
                if (yvec[i] > ymax1) ymax1 = yvec[i];
                if (yvec[i] < ymin1) ymin1 = yvec[i];
            }
            xdenom = xmax1 - xmin1;
            ydenom = ymax1 - ymin1;
            if (xdenom == 0 || ydenom == 0) {
                System.out.println("*** ERROR PlotContoursRect(): "+
                                   "bad input data: max=min for x or y.");
                System.exit(0);
            }
        } else {
            xmin1=xmin;
            xmax1=xmax;
            ymin1=ymin;
            ymax1=ymax;
        }
        
        for (i = 0; i < nx; i++)
            xvec[i] = (xvec[i] - xmin1) / xdenom;
        for (i = 0; i < ny; i++)
            yvec[i] = (yvec[i] - ymin1) / ydenom;

        for (i=0; i<ncontours; i++) {
            zvalue = contours[i];
            ContourReturn c = get_contour_line(z,xvec,yvec,nx,ny,zvalue);
            x = c.getX();
            y = c.getY();
            n = c.n;
            // Need to break up the line into multiple 500-pt lines:
            //  ... is n-2 b/c last pt in 500 line blocks are done in next
            //     block too; 
            //   i.e. 1002 pts need 3 blocks, but 1001 only needs 2 blocks.

            for (j=0; j<n; j+=2) {
                g.drawLine((int)(x[j]*xscale),(int)(y[j]*yscale),
                           (int)(x[j+1]*xscale),(int)(y[j+1]*yscale));
            }
        }

        //  Now unscale the input vectors:
        for (i = 0; i < nx; i++)
            xvec[i] = xdenom * xvec[i] + xmin1;
        for (i = 0; i < ny; i++)
            yvec[i] = ydenom * yvec[i] + ymin1;
        return;
    }


    // Find all line segments in a contour line.
    // USES:
    //     get_segments
    // CALLED FROM:
    //     plot_contours_polar, plot_contours_rect
    
    // Returns x,y coords of iso-z-line specified by the given zvalue.
    // Also returns number of points in this line: n.
    /*  The line is actually many segments: moveto/lineto.
        Since there could be more than one "connected line segment"
        for any zvalue, we just return these segments.
        Figuring out how the segments are arranged into
        connected line segments is not attempted.  */

    private ContourReturn get_contour_line(double z[], 
                                           double xvec[], double yvec[],
                                           int nx, int ny, double zvalue) {
        int i, j, k;
        int nmax;
        double xs[] = new double[10];
        double ys[] = new double[10];
        int nsegs;
        double x1, x2, y1, y2;
        double z11, z12, z21, z22;

        // OUTPUT PARAMETERS,  INSIDE ContourReturn
        //double x[], y[];  // actually Vectors
        //int n;
        ContourReturn c = new ContourReturn();
        
        nmax = NMAX;
        c.n = -1;   // was 0

        for (i=0; i<(nx-1); i++) {
            x1 = xvec[i];
            x2 = xvec[i+1];
            for (j=0; j<(ny-1); j++) {
                y1 = yvec[j];
                y2 = yvec[j+1];
                z11 = z[(j-1)*nx + i];
                z12 = z[(j)*nx + i];
                z21 = z[(j-1)*nx + i+1];
                z22 = z[(j)*nx + i+1];

                //get_segments(x1,x2,y1,y2,z11,z12,z21,z22,zvalue,
                //             xs,ys,nsegs);
                SegmentReturn s = get_segments(x1,x2,y1,y2,
                                               z11,z12,z21,z22,zvalue);

                for (k=0; k<s.nsegs; k+=2) {
                    c.n++;
                    if (c.n < nmax) {
                        //x[n] = xs[k];
                        //y[n] = ys[k];
                        c.x.addElement(new Double(s.xs[k]));
                       c.y.addElement(new Double(s.ys[k]));
                    } else return c;
                    c.n++;
                    if (c.n < nmax) {
                        //x[n] = xs[k+1];
                        //y[n] = ys[k+1];
                        c.x.addElement(new Double(s.xs[k+1]));
                        c.y.addElement(new Double(s.ys[k+1]));
                    } else return c;
                }
            }
        }
        return c;
    }


    // Find all line segments in a pixel.
    //
    // CALLED FROM:
    //     get_contour_line
    //
    /* ! A pixel has 4 edges, defined by the 4 variables {\tt x1}, {\tt x2},
       {\tt y1}, and {\tt y2}.
       The values of the function at each corner of the pixel are defined
       by the variables {\tt z11}, {\tt z12}, {\tt z21}, and {\tt z22}.
       The first subscript refers to the x-coordinate, the second to the y-coordinate.
       This is shown in the figure below:

                    z12            z22
       y2         *---------------*
                  |               |
                  |               |
                  |               |
                  |               |
                  |               |
                  | z11           | z21
       y1         *---------------*
       
                  x1              x2
       
       If the z-value is between any pair of z's on the edges, we note where
       along that edge that would occur.
       If there are 2 of these, we connect them with a line segment and return it.
       If there are 3 or 4 of these, we connect them, somewhat arbitrarily,
       and return multiple line segments.
    */
    
    private SegmentReturn get_segments(double x1,double x2,double y1,double y2,
                                       double z11,double z12,double z21,double z22,
                                       double zvalue) {

        // OUTPUT STUFF:  inside SegmentReturn
        //double xs[] = new double[10];
        //double ys[] = new double[10];
        //int nsegs;
        // END OUTPUT STUFF
        SegmentReturn s = new SegmentReturn();

        double zmin, zmax, m1, m2;
        int n_intercepts;
        double xint[] = new double[10];
        double yint[] = new double[10];
        double frac;

        m1 = Math.min(z11, z12);
        m2 = Math.min(z21, z22);
        zmin = Math.min(m1, m2);
        m1 = Math.max(z11, z12);
        m2 = Math.max(z21, z22);
        zmax = Math.max(m1, m2);

        if (zvalue < zmin || zvalue > zmax) {
            // contour not possible
            s.nsegs = 0;
            return s;
        }

        // go through each edge of rectangle, finding intercepts with zvalue:
        n_intercepts=0;

        // 1. bottom edge:
        zmin = Math.min(z11, z21);
        zmax = Math.max(z11, z21);
        if (zvalue >= zmin && zvalue <= zmax) {
            n_intercepts= n_intercepts+1;
            frac=(x2-x1)*(zvalue-zmin)/(zmax-zmin);
            if (Math.abs(zmin-z11) <= 0.1*Math.abs(zmin)) {
                // zmin is at x1:
                xint[n_intercepts] = x1 + frac;
            } else {
                // zmin is at x2:
                xint[n_intercepts] = x2 - frac;
            }
            yint[n_intercepts] = y1;
        }

        // 2. left edge:
        zmin = Math.min(z11,z12);
        zmax = Math.max(z11,z12);
        if (zvalue >= zmin && zvalue <= zmax) {
            n_intercepts= n_intercepts+1;
            frac=(y2-y1)*(zvalue-zmin)/(zmax-zmin);
            if (Math.abs(zmin-z11) <= 0.1*Math.abs(zmin)) {
                // zmin is at y1:
                yint[n_intercepts] = y1 + frac;
            } else {
                // zmin is at y2:
                yint[n_intercepts] = y2 - frac;
            }
            xint[n_intercepts] = x1;
        }

        // 3. top edge:
        zmin = Math.min(z12,z22);
        zmax = Math.max(z12,z22);
        if (zvalue >= zmin && zvalue <= zmax) {
            n_intercepts= n_intercepts+1;
            frac=(x2-x1)*(zvalue-zmin)/(zmax-zmin);
            if (Math.abs(zmin-z12) <= 0.1*Math.abs(zmin)) {
                // zmin is at x1:
                xint[n_intercepts] = x1 + frac;
            } else {
                // zmin is at x2:
                xint[n_intercepts] = x2 - frac;
            }
            yint[n_intercepts] = y2;
        }

        // 4. right edge:
        zmin = Math.min(z21,z22);
        zmax = Math.max(z21,z22);
        if (zvalue >= zmin && zvalue <= zmax) {
            n_intercepts= n_intercepts+1;
            frac=(y2-y1)*(zvalue-zmin)/(zmax-zmin);
            if (Math.abs(zmin-z21) <= 0.1*Math.abs(zmin)) {
                // zmin is at y1:
                yint[n_intercepts] = y1 + frac;
            } else {
                // zmin is at y2:
                yint[n_intercepts] = y2 - frac;
            }
            xint[n_intercepts] = x2;
        }

        if (n_intercepts == 0) {
            s.nsegs = 0;
        } else if (n_intercepts == 1) {
            // NOT POSSIBLE
            s.nsegs = 0;
        } else if (n_intercepts == 2) {
            // connect them with a line:
            s.nsegs = 1;
            s.xs[0] = xint[0];
            s.xs[1] = xint[1];
            s.ys[0] = yint[0];
            s.ys[1] = yint[1];
        } else if (n_intercepts == 3) {
            // connect as 2 lines (a guess)
            s.nsegs = 2;
            s.xs[0] = xint[0];
            s.xs[1] = xint[1];
            s.ys[0] = yint[0];
            s.ys[1] = yint[1];

            s.xs[2] = xint[1];
            s.xs[3] = xint[2];
            s.ys[2] = yint[1];
            s.ys[3] = yint[2];
        }else if (n_intercepts == 4) {
            // connect as 3 lines (a guess)
            s.nsegs = 3;
            s.xs[0] = xint[0];
            s.xs[1] = xint[1];
            s.ys[0] = yint[0];
            s.ys[1] = yint[1];

            s.xs[2] = xint[1];
            s.xs[3] = xint[2];
            s.ys[2] = yint[1];
            s.ys[3] = yint[2];

            s.xs[4] = xint[2];
            s.xs[5] = xint[3];
            s.ys[4] = yint[2];
            s.ys[5] = yint[3];
        }
        return s;
    }

}


class ContourReturn {
    public Vector<Double> x;
    public Vector<Double> y;
    int n;

    public ContourReturn() {
        super();
        this.x = new Vector<Double>();
        this.y = new Vector<Double>();
        this.n = 0;
    }

    
    public double[] getX() {
        double newX[];
        n = x.size();
        newX = new double[n];
        for (int i=0; i<n; i++)
            newX[i] = ((Double)x.elementAt(i)).doubleValue();
        return newX;
    }
    
    public double[] getY() {
        double newY[];
        n = y.size();
        newY = new double[n];
        for (int i=0; i<n; i++)
            newY[i] = ((Double)y.elementAt(i)).doubleValue();
        return newY;
    }


}


class SegmentReturn{
    double xs[] = new double[10];
    double ys[] = new double[10];
    int nsegs;

    public SegmentReturn() {
        super();
    }
}
