// toro_area.java   solution to PDE

import java.awt.*;
import java.awt.event.*;

public class toro_area extends Frame
{
  int height = 600;
  int width = 600;
  // specific plot boundaries
  double x, y, z, t;         // orthographic z up, x across, y back
  double xang = Math.PI/4.0; // 45 degrees
  double yang = Math.PI/6.0; // 30 degrees
  double r1 = 2.0;
  double r2 = 0.5;
  double xmin = -(r1+r2);
  double xmax = r1+r2;
  double ymin = -(r1+r2);
  double ymax = r1+r2;
  double zmin = -(r1+r2); // actually -r2, keep equal scaling in plot
  double zmax = (r1+r2);  // actually  r2
  // surface generation
  int nth = 96;
  int nph = 24;
  int n = nth*nph; // total points
  double far[] = new double[n];
  double ftheta[] = new double[n];
  double fphi[] = new double[n];
  int ixk[] = new int[n];

  double dtheta = 2.0*Math.PI/(double) nth;
  double dphi = 2.0*Math.PI/(double)nph;
  double theta_max = 2.0*Math.PI;
  double phi_max = 2.0*Math.PI;
  double area, tarea, rerr;
  double volume, tvolume, verr;

  toro_area()
  {
    double x1, x2, x3, x4;
    double y1, y2, y3, y4;
    double z1, z2, z3, z4;
    double theta, phi;
    int k;
    heapsort H = new heapsort();

    setTitle("toro_area");
    setSize(width,height);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
    // build points, sort for display
    for(int i=0; i<nth; i++)
    {
      theta = dtheta*(double)i;
      for(int j=0; j<nph; j++)
      {
	phi = dphi*(double)j;
	x1 = fx(r1,r2,phi,theta);
	x2 = fx(r1,r2,phi,theta+dtheta);
        x3 = fx(r1,r2,phi+dphi,theta+dtheta);
        x4 = fx(r1,r2,phi+dphi,theta);
	y1 = fy(r1,r2,phi,theta);
	y2 = fy(r1,r2,phi,theta+dtheta);
        y3 = fy(r1,r2,phi+dphi,theta+dtheta);
        y4 = fy(r1,r2,phi+dphi,theta);
	z1 = fz(r1,r2,phi,theta);
	z2 = fz(r1,r2,phi,theta+dtheta);
        z3 = fz(r1,r2,phi+dphi,theta+dtheta);
        z4 = fz(r1,r2,phi+dphi,theta);
        k = i*nph+j;
        x1 = (x1+x2+x3+x4)/4.0;
        y1 = (y1+y2+y3+y4)/4.0;
        far[k] = y1*Math.sin(yang)-x1*Math.sin(xang);
        far[k] = -far[k];
        ftheta[k] = theta;
        fphi[k] = phi;
        // System.out.println("k="+k+", far="+far[k]+", th="+ftheta[k]+
        //                    ", ph="+fphi[k]);
      } // phi
    } // theta
    H.sortdindex(far, n, ixk);
    // System.out.println(ixk[0]+" "+ixk[1]+" "+ixk[2]+" "+
    //                    ixk[n-2]+" "+ixk[n-1]);

  }

  public void paint(Graphics g)
  {
    double x1, x2, x3, x4;
    double y1, y2, y3, y4;
    double z1, z2, z3, z4;
    int xi, yi, xj, yj;
    int xpts[] = new int[4];
    int ypts[] = new int[4];
    double phi, theta;
    double l12, l23, l34, l41;

    System.out.println("toro_area.java  r1="+r1+", r2="+r2+
                       ", ntheta="+nth+", nphi="+nph);
    g.setColor(Color.black);
    g.drawString("toro_area.java  r1="+r1+", r2="+r2+
                 ", ntheta="+nth+", nphi="+nph, 20, 50);
    g.drawString("0<theta<2Pi sweeps r1, 0<phi<2Pi sweeps r2", 20, 70);
    g.drawString("x = (r1+r2*sin(phi))*cos(theta)", 20, 90);
    g.drawString("y = (r1+r2*sin(phi))*sin(theta)", 20, 110);
    g.drawString("z = (r2*cos(phi)", 20, 130);
    g.drawString("plotted orthogonal   Xangle="+xang+", Yangle="+yang,
                 20, 160);
    g.drawString("Xplot = x*cos(Xangle)+y*cos(Yangle)    then scaled",
                 20, 180);
    g.drawString("Yplot = x*sin(Xangle)-y*sin(Yangle)+z  then scaled",
                 20, 200);

    area = 0.0;
    volume = 0.0;
    for(int k=0; k<n; k++)
    {
      theta = ftheta[ixk[k]];
      phi = fphi[ixk[k]];
      x1 = fx(r1,r2,phi,theta);
      x2 = fx(r1,r2,phi,theta+dtheta);
      x3 = fx(r1,r2,phi+dphi,theta+dtheta);
      x4 = fx(r1,r2,phi+dphi,theta);
      y1 = fy(r1,r2,phi,theta);
      y2 = fy(r1,r2,phi,theta+dtheta);
      y3 = fy(r1,r2,phi+dphi,theta+dtheta);
      y4 = fy(r1,r2,phi+dphi,theta);
      z1 = fz(r1,r2,phi,theta);
      z2 = fz(r1,r2,phi,theta+dtheta);
      z3 = fz(r1,r2,phi+dphi,theta+dtheta);
      z4 = fz(r1,r2,phi+dphi,theta);
      // normal sets color
      normal(x1,y1,z1, x2,y2,z2, x3,y3,z3, g);
      xpts[0] = xortho(x1, y1, z1);
      ypts[0] = yortho(x1, y1, z1);
      xpts[1] = xortho(x2, y2, z2);
      ypts[1] = yortho(x2, y2, z2);
      xpts[2] = xortho(x3, y3, z3);
      ypts[2] = yortho(x3, y3, z3);
      xpts[3] = xortho(x4, y4, z4);
      ypts[3] = yortho(x4, y4, z4);
      g.fillPolygon(xpts, ypts, 4);

      l12 = side(x1, y1, z1, x2, y2, z2);
      l23 = side(x2, y2, z2, x3, y3, z3);
      l34 = side(x3, y3, z3, x4, y4, z4);
      l41 = side(x4, y4, z4, x1, y1, z1);
      // if(theta+dtheta<=theta_max && phi+dphi<=phi_max)
      area = area + ((l12+l34)/2.0)*((l23+l41)/2.0); // approx
      volume = volume + ((l12+l34)/2.0)*((l23+l41)/2.0)*r2/2.0;
                                 // approximation 
    } // end k
    System.out.println("computed area="+area);
    System.out.println("computed volume="+volume);
    tarea = 2.0*Math.PI*r1 * 2.0*Math.PI*r2;
    tvolume = 2.0*Math.PI*r1 * Math.PI*r2*r2;
    rerr = (tarea-area)/tarea;
    verr = (tvolume-volume)/tvolume;
    System.out.println("theory area="+tarea+", relative error="+rerr);
    g.setColor(Color.black);
    g.drawString("colors based on normal to surface: x red, y green, z blue",
                 20, 510);
    g.drawString("computed area="+area, 20, 530);
    g.drawString("theory area="+tarea+", relative error="+rerr, 20, 550);
    g.drawString("computed volume="+volume, 20, 570);
    g.drawString("theory volume="+tvolume+", relative error="+verr, 20, 590);

  } // end paint

  double fx(double ra, double rb, double phi, double theta) // function x to be plotted
  {
      return (ra+rb*Math.sin(phi))*Math.cos(theta);
  }

  double fy(double ra, double rb, double phi, double theta) // function y to be plotted
  {
      return (ra+rb*Math.sin(phi))*Math.sin(theta);
  }


  double fz(double ra, double rb, double phi, double theta) // function z to be plotted
  {
      return rb*Math.cos(phi);
  }

  int xortho(double x1, double y1, double z1) // screen X projection
  {
      double xp, yp, zp;
      xp = (x1-xmin)/(xmax-xmin);
      yp = (y1-ymin)/(ymax-ymin);
      zp = (z1-zmin)/(zmax-zmin);
      return 30+(int)(width*0.57*(xp*Math.cos(xang)+yp*Math.cos(yang)));
  } // end xortho

  int yortho(double x1, double y1, double z1) // screen Y projection
  {
      double xp, yp, zp;
      xp = (x1-xmin)/(xmax-xmin);
      yp = (y1-ymin)/(ymax-ymin);
      zp = (z1-zmin)/(zmax-zmin);
      return height-150-(int)(height*0.38*(-xp*Math.sin(xang)+
                                            yp*Math.sin(yang)+zp));
  } // end yortho

  void normal(double x0, double y0, double z0,
              double x1, double y1, double z1,
              double x2, double y2, double z2,
              Graphics g)
  {
    double ax, ay, az, bx, by, bz, nx, ny, nz, denom;
    float red, grn, blu; 
    ax = x2 - x0;
    ay = y2 - y0;
    az = z2 - z0;
    bx = x1 - x0;
    by = y1 - y0;
    bz = z1 - z0;
    nx = ay*bz-az*by;
    ny = az*bx-ax*bz;
    nz = ax*by-ay*bx;

    denom = Math.sqrt(nx*nx+ny*ny+nz*nz);
    nx = nx/denom;
    ny = ny/denom;
    nz = nz/denom;
    // System.out.println(x0+y0+z0);
    // System.out.println(x1+y1+z1);
    // System.out.println(x2+y2+z2);
    // System.out.println("nx="+nx+", ny="+ny+", nz="+nz);
    red = (float)((nx+1.0)/2.0); // change color based on normal vector
    grn = (float)((ny+1.0)/2.0); // normal from -1 to +1, color 0 to 1
    blu = (float)((nz+1.0)/2.0); // must be at least 2
    g.setColor(new Color(red,grn,blu));
  } // end normal

  double side(double x1, double y1, double z1, double x2, double y2, double z2)
  {
    // length of one side
    return Math.sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
  } // end side

  public static void main(String args[])
  {
    new toro_area();
  }
} // end toro_area.java