// draw_sphere_deriv.java  sphereical coordinates
// x = r cos(theta) sin(phi)    0 < theta < 2 Pi
// y = r sin(theta) sin(phi)    0 < phi < Pi restricted for derivatives
// z = r cos(phi)               r < 0        restricted for derivatives
// r = sqrt(x^2+y^2+z^2)
// theta = arctan(y,x) if theta<0 then theta = 2 Pi - theta
// phi = arctan(sqrt(x^2+y^2),z)
//
// compute del U(radius,theta,phi) r=radius, t=theta angle, p=phi angle
// del U  =  (d/dr U, 1/(r*sin(p)) d/dt U, (1/r) d/dp U)
//
// compute del^2 U(radius,theta,phi) r=radius, t=theta angle, p=phi angle
// del^2 U  = (1/r^2) d/dr(r^2 d/dr U) + 1/(r^2 sin(p)^2) d^2/dt^2 U +
//            (1/(r^2 sin(p) d/dp(sin(p) d/dp U)
//          = (2/r) d/dr U + d^2/dr^2 U + (1/(r^2 sin(p)^2) d^2/dt^2 U +
//            (cos(p)/(r^2 sin(t)) d/dp U + (1/(r^2)) d^2/dp^2 U

import java.awt.*;
import java.awt.event.*;
import java.text.*;
import java.io.*;

// uses nuderiv,java
public class draw_sphere_deriv extends Frame
{
  double Pi = 3.1415926535897932384626433832795028841971;
  int Nr = 5;
  double Rmin = 0.3;
  double Rmax = 1.0;
  int Nt = 7;
  double Tmin = 0.1;
  double Tmax = 300.0*Pi/180.0;
  int Np = 5;
  double Pmin = 0.2;
  double Pmax = Pi-0.5;
  double Rg[] = new double[Nr]; // r grid values
  double Tg[] = new double[Nt]; // theta grid values
  double Pg[] = new double[Np]; // p grid values
  int Nrtp = Nr*Nt*Np;
  double Ua[] = new double[Nrtp];
  double Cr[] = new double[Nr]; // first deriv wrt r
  double Crr[] = new double[Nr]; // second deriv wrt r
  double Ct[] = new double[Nt]; // first deriv wrt t
  double Ctt[] = new double[Nt]; // second deriv wrt t
  double Cp[] = new double[Np]; // first deriv wrt p
  double Cpp[] = new double[Np]; // second deriv wrt p
  double R, T, P, Ana, Cmp, Err;
  double Urv, Urrv, Utv, Uttv, Upv, Uppv; // computed derivatives
  double x, y, z, rtr, ttr, ptr;
  double Maxerrur=0.0, Maxerrurr=0.0, Maxerrut=0.0, Maxerrutt=0.0;
  double Maxerrup=0.0, Maxerrupp=0.0, Maxerr=0.0, Maxerrlap=0.0;

  // for plotting
  int width = 500;
  int height = 500;
  double xmin = -1.0;
  double xmax =  1.0;
  double ymin = -1.0;
  double ymax =  1.0;
  double zmin = -1.0;
  double zmax =  1.0;

  double ya = 0.5;     //  orthographic typ 0.1 to 0.5 
  double xs = (0.666*(width-100))/(xmax-xmin);
  double ys = (0.666*(height-100))/(ymax-ymin);
  int xp, yp; // computed by parab

  draw_sphere_deriv()
  {
    System.out.println("draw_sphere_deriv.java starting");
    System.out.println("r gird values, not equally spaced");
    Rg[0] = 0.3;
    Rg[1] = 0.5;
    Rg[2] = 0.7;
    Rg[3] = 0.9;
    Rg[4] = 1.0;
    for(int I=0; I<Nr; I++)
       System.out.println("Rg["+(I)+"]="+(Rg[I]));
    System.out.println("theta gird values, not equally spaced");
    Tg[0] = 0.1;
    Tg[1] = 40.0*Pi/180.0;
    Tg[2] = 85.0*Pi/180.0;
    Tg[3] = 130.0*Pi/180.0;
    Tg[4] = 185.0*Pi/180.0;
    Tg[5] = 240.0*Pi/180.0;
    Tg[6] = 300.0*Pi/180.0;
    for(int J=0; J<Nt; J++)
       System.out.println("Tg["+(J)+"]="+(Tg[J]));
    System.out.println("phi gird values, not equally spaced");
    Pg[0] = 0.2;
    Pg[1] = 0.9;
    Pg[2] = 1.7;
    Pg[3] = 2.4;
    Pg[4] = Pi-0.5;
    for(int K=0; K<Np; K++)
       System.out.println("Pg["+(K)+"]="+(Pg[K]));

    System.out.println("analytic partial second derivative");
    for(int I=0; I<Nr; I++)
      for(int J=0; J<Nt; J++)
	for(int K=0; K<Np; K++)
	{
          Ua[S(I,J,K)] = U(Rg[I],Tg[J],Pg[K]);
          System.out.println("Ua["+(I)+","+(J)+","+(K)+
                             "]="+(Ua[S(I,J,K)])+
                             ", D2="+(Uuu(Rg[I],Tg[J],Pg[K])));
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          if(I==0 && J==0 && K==0) 
            System.out.println("x="+x+", y="+y+", z="+z);
          rtr = Math.sqrt(x*x+y*y+z*z);
          ttr = Math.atan2(y, x);
          if(ttr<0.0) ttr = 2.0*Pi + ttr; // atan2 -Pi to Pi fix
          ptr = Math.atan2(Math.sqrt(x*x+y*y),z);
          Err = Math.abs(rtr-Rg[I]);
          if(Err>0.001) 
            System.out.println("R="+Rg[I]+", rtr="+rtr+", err="+Err);
          Err = Math.abs(ttr-Tg[J]);
          if(Err>0.001) 
	    {System.out.println("T="+Tg[J]+", rtr="+ttr+", err="+Err);
	     System.out.println("x="+x+", y="+y+",z="+z);}
          Err = Math.abs(ptr-Pg[K]);
          if(Err>0.001) 
            System.out.println("P="+Pg[K]+", ptr="+ptr+", err="+Err);
          Err = Math.abs(U(Rg[I],Tg[J],Pg[K])-Uxyz(x,y,z));
          if(Err>0.001) 
	    {System.out.println("U="+U(Rg[I],Tg[J],Pg[K])+
                                ", Uxyz="+Uxyz(x,y,z)+", err="+Err);
	     System.out.println("x="+x+", y="+y+",z="+z);}

	}

    for(int I=0; I<Nr; I++)
    {
      R = Rg[I];
      for(int J=0; J<Nt; J++)
      {
        T = Tg[J];
        for(int K=0; K<Np; K++)
	{
          P = Pg[K];
          // compute Ur, Urr, Utt, Upp and call D2
          new nuderiv(1, Nr, I, Rg, Cr);
          new nuderiv(2, Nr, I, Rg, Crr);
          Urv = 0.0;
          Urrv = 0.0;
          for(int Ii=0; Ii<Nr; Ii++)
	  {
            Urv = Urv + Cr[Ii]*U(Rg[Ii],T,P);
            Urrv = Urrv + Crr[Ii]*U(Rg[Ii],T,P);
          }
          Utv = 0.0;
          Uttv = 0.0;
          new nuderiv(1, Nt, J, Tg, Ct);
          new nuderiv(2, Nt, J, Tg, Ctt);
          for(int Jj=0; Jj<Nt; Jj++)
	  {
            Utv = Utv + Ct[Jj]*U(R,Tg[Jj],P);
            Uttv = Uttv + Ctt[Jj]*U(R,Tg[Jj],P);
          }
          Upv = 0.0;
          Uppv = 0.0;
          new nuderiv(1, Np, K, Pg, Cp);
          new nuderiv(2, Np, K, Pg, Cpp);
          for(int Kk=0; Kk<Np; Kk++)
	  {
            Upv = Upv + Cp[Kk]*U(R,T,Pg[Kk]);
            Uppv = Uppv + Cpp[Kk]*U(R,T,Pg[Kk]);
          }

          Cmp = D2(Urv, Urrv, Uttv, Upv, Uppv, R);
          Ana = Uuu(R,T,P);
          Err = Ana - Cmp;
          Maxerr = Math.max(Maxerr, Err);
          System.out.println("i="+(I)+
                  ", j="+(J)+
                  ", k="+(K)+
                  ", Ana="+(Ana)+
                  ", Cmp="+(Cmp)+
                  ", err="+(Err));
          Err = Ur(R,T,P)-Urv;
          Maxerrur = Math.max(Maxerrur, Err);
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Ur="+(Ur(R,T,P))+
                  ", Urv="+(Urv)+
                  ", err="+(Err));
          Err = Urr(R,T,P)-Urrv;
          Maxerrurr = Math.max(Maxerrurr, Err);
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Urr="+(Urr(R,T,P))+
                  ", Urrv="+(Urrv)+
                  ", err="+(Err));
          Err = Ut(R,T,P)-Utv;
          Maxerrut = Math.max(Maxerrut, Err);
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Ut="+(Ut(R,T,P))+
                  ", Utv="+(Utv)+
                  ", err="+(Err));
          Err = Utt(R,T,P)-Uttv;
          Maxerrutt = Math.max(Maxerrutt, Err);
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Utt="+(Utt(R,T,P))+
                  ", Uttv="+(Uttv)+
                  ", err="+(Err));
          Err = Up(R,T,P)-Upv;
          Maxerrupp = Math.max(Maxerrupp, Err);
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Up="+(Upp(R,T,P))+
                  ", Upv="+(Upv)+
                  ", err="+(Err));
          System.out.println(" ");
          Err = Upp(R,T,P)-Uppv;
          Maxerrupp = Math.max(Maxerrupp, Err);
          System.out.println("i="+(I)+", j="+(J)+", k="+(K)+
                  ", Upp="+(Upp(R,T,P))+
                  ", Uppv="+(Uppv)+
                  ", err="+(Err));
          System.out.println(" ");
        }
      }
    }
    System.out.println("Maxerr Ur="+(Maxerrur));
    System.out.println("Maxerr Urr="+(Maxerrurr));
    System.out.println("Maxerr Ut="+(Maxerrut));
    System.out.println("Maxerr Utt="+(Maxerrutt));
    System.out.println("Maxerr Up="+(Maxerrup));
    System.out.println("Maxerr Upp="+(Maxerrupp));
    System.out.println("Maxerr deriv="+(Maxerr));

    write_file();

    setTitle("draw_sphere_deriv");
    setSize(width,height);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
    repaint();

    System.out.println("draw_sphere_deriv.java finishing");
  } //end draw_sphere_deriv

  double Uxyz(double x, double y, double z)
  {
    double rtr, ttr, ptr;
    rtr = Math.sqrt(x*x+y*y+z*z);
    ttr = Math.atan2(y, x);
    if(ttr<0.0) ttr = 2.0*Pi + ttr; // atan2 -Pi to Pi fix
    ptr = Math.atan2(Math.sqrt(x*x+y*y),z);
    return rtr*rtr*rtr + ttr*ttr*ttr + ptr*ptr*ptr;
  } // end Uxyz

  double U(double R, double T, double P)
  {
    return R*R*R + T*T*T + P*P*P;
  } // end U

  double Ur(double R, double T, double P)
  {
    return 3.0*R*R;
  } // end Ur

  double Urr(double R, double T, double P)
  {
    return 6.0*R;
  } // end Urr;

  double Ut(double R, double T, double P)
  {
    return 3.0*T*T;
  } // end Ut

  double Utt(double R, double T, double P)
  {
    return 6.0*T;
  } // end Utt;

  double Up(double R, double T, double P)
  {
    return 3.0*P*P;
  } // end Up

  double Upp(double R, double T, double P)
  {
    return 6.0*P;
  } // end Upp

  double Uuu(double R, double T, double P)
  {
    return 2.0*Ur(R,T,P)/R + Urr(R,T,P) + 
	   Utt(R,T,P)/(R*R*Math.sin(P)*Math.sin(P)) +
           Math.cos(P)*Up(R,T,P)/(R*R*Math.sin(P)) + 
	   Upp(R,T,P)/(R*R);
  } // end Uuu

  double D2(double Urv, double Urrv, double Uttv, double Upv,
            double Uppv, double R)
  {
    return 2.0*Urv/R + Urrv + 
	   Uttv/(R*R*Math.sin(P)*Math.sin(P)) + 
	   Math.cos(P)*Upv/(R*R*Math.sin(P)) + 
           Uppv/(R*R);
  } // end D2

  int S(int I, int J, int K)
  { // index of each grid point
    return I*Nt*Np + J*Np + K;
  } // end S

  // plotting functions

  public void paint(Graphics g)
  {
    double x, y, z;
    double z1, z2, hz;
    int xp1, yp1, xp2, yp2; // line segment

    // labels
    set_color(0, g);
    g.drawString("draw_sphere_deriv", 200, 50);
    g.drawString("xmin", 25, height-15);
    g.drawString("xmax", (int)(0.66*width-70), height-15);
    g.drawString("ymin", (int)(0.66*width-30), height-25);
    g.drawString("ymax", (int)(width-80), (int)(0.66*height+25));
    g.drawString("zmin", (int)(width-75), (int)(0.66*height+5));
    g.drawString("zmax", (int)(width-75), (int)(85));

    // draw box, bottom, top, sides
    set_color(0, g);
    ortho(xmin, ymin, zmin); xp1=xp; yp1=yp;
    ortho(xmax, ymin, zmin); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, zmin); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymax, zmin); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymin, zmin); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymin, zmax); xp1=xp; yp1=yp;
    ortho(xmax, ymin, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, zmax); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymax, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymin, zmax); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);

    ortho(xmin, ymin, zmin); xp1=xp; yp1=yp;
    ortho(xmin, ymin, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymax, zmin); xp1=xp; yp1=yp;
    ortho(xmin, ymax, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymin, zmin); xp1=xp; yp1=yp;
    ortho(xmax, ymin, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, zmin); xp1=xp; yp1=yp;
    ortho(xmax, ymax, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);

    // plot function Pg
    set_color(2, g);
    g.drawString("Phi", 25, 40);
    for(int I=0; I<Nr; I++)
      for(int J=0; J<Nt; J++)
	for(int K=0; K<Np-1; K++)
	{
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          ortho(x, y, z); xp1=xp; yp1=yp;
	  if(K+1<Np)
	  {
	    x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K+1]);
	    y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K+1]);
	    z = Rg[I]*Math.cos(Pg[K+1]);
            ortho(x, y, z); xp2=xp; yp2=yp;
	  }
	  else
	  {
	    break;
	  }
          g.drawLine(xp1, yp1, xp2, yp2);
	}

    // plot function Rg
    set_color(5, g);
    g.drawString("Radius", 25, 60);
    for(int J=0; J<Nt; J++)
      for(int K=0; K<Np; K++)
        for(int I=0; I<Nr-1; I++)
	{
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          ortho(x, y, z); xp1=xp; yp1=yp;
	  if(I+1<Nr)
	  {
	    x = Rg[I+1]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	    y = Rg[I+1]*Math.sin(Tg[J])*Math.sin(Pg[K]);
	    z = Rg[I+1]*Math.cos(Pg[K]);
	  }
	  else
	  {
	    x = Rg[0]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	    y = Rg[0]*Math.sin(Tg[J])*Math.sin(Pg[K]);
	    z = Rg[0]*Math.cos(Pg[K]);
	  }
          ortho(x, y, z); xp2=xp; yp2=yp;
          g.drawLine(xp1, yp1, xp2, yp2);
	}


    // plot function Tg
    set_color(6, g);
    g.drawString("Theta", 25, 80);
    for(int K=0; K<Np; K++)
      for(int I=0; I<Nr; I++)
        for(int J=0; J<Nt-1; J++)
	{
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          ortho(x, y, z); xp1=xp; yp1=yp;
	  if(J+1<Nt)
	  {
	    x = Rg[I]*Math.cos(Tg[J+1])*Math.sin(Pg[K]);
	    y = Rg[I]*Math.sin(Tg[J+1])*Math.sin(Pg[K]);
	    z = Rg[I]*Math.cos(Pg[K]);
	  }
	  else
	  {
	    x = Rg[I]*Math.cos(Tg[0])*Math.sin(Pg[K]);
	    y = Rg[I]*Math.sin(Tg[0])*Math.sin(Pg[K]);
	    z = Rg[I]*Math.cos(Pg[K]);
	  }
          ortho(x, y, z); xp2=xp; yp2=yp;
          g.drawLine(xp1, yp1, xp2, yp2);
	}

  } // end paint

  // orthographic projection, sets xp,yp
  void ortho(double x, double y, double z)
  {
    double x1, y1, z1, x2, y2, x3, y3;
    x1 = x-xmin;   // zero base
    y1 = y-ymin; 
    z1 = z-zmin;
    x2 = x1+y1*ya; // ortho z is up, y is back
    y2 = z1+y1*ya;
    x3 = x2*xs;    // scale 
    y3 = y2*ys;
    x3 = x3+25;
    y3 = height-y3-25; // away from boarder
    xp = (int)x3;
    yp = (int)y3;
  }

  void set_color(int i, Graphics g)
  {
    if(i<0) {g.setColor(Color.black); return;}
    if(i>9) {g.setColor(Color.gray); return; }
    switch(i)
    {
      case 0: g.setColor(Color.black); break;
      case 1: g.setColor(Color.pink); break; // should be brown
      case 2: g.setColor(Color.red); break;
      case 3: g.setColor(Color.orange); break;
      case 4: g.setColor(Color.yellow); break;
      case 5: g.setColor(Color.green); break;
      case 6: g.setColor(Color.blue); break;
      case 7: g.setColor(Color.magenta); break;
      case 8: g.setColor(Color.cyan); break;
      case 9: g.setColor(Color.gray); break;
    }
  } // end set_color

  void write_file()
  {
    BufferedWriter outp;
    PrintWriter file_out;
    double x, y, z;
    double rgv, tgv, pgv; // steping values for plot

    try
    {
      outp = new BufferedWriter(new FileWriter("spherical_deriv.plot"));
      file_out = new PrintWriter(outp);
    }
    catch(Exception e)
    {
      System.out.println("could not open spherical_deriv.plot for write ");
      return;
    }
    file_out.println("set hidden3d ");
    file_out.println("set xrange ["+xmin+":"+xmax+"]");
    file_out.println("set yrange ["+ymin+":"+ymax+"]");
    file_out.println("splot \"spherical_deriv.dat\" with lines ");
    file_out.println("exit");
    file_out.close();

    try
    {
      outp = new BufferedWriter(new FileWriter("spherical_deriv.dat"));
      file_out = new PrintWriter(outp);
    }
    catch(Exception e)
    {
      System.out.println("can't open spherical.dat for write");
      return;
    }
    for(int I=0; I<Nr; I++)
    {
      for(int J=0; J<Nt; J++)
      {
        for(int K=0; K<Np-1; K++)
        {
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          file_out.println(x+" "+y+" "+z);
          pgv = Pg[K];
	  if(K+1<Np)
	  {
	    while(Pg[K+1]>pgv+0.1)
	    {
	      pgv = pgv + 0.1;
	      x = Rg[I]*Math.cos(Tg[J])*Math.sin(pgv);
	      y = Rg[I]*Math.sin(Tg[J])*Math.sin(pgv);
	      z = Rg[I]*Math.cos(pgv);
              file_out.println(x+" "+y+" "+z);
	    } // end while
	    x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K+1]);
	    y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K+1]);
	    z = Rg[I]*Math.cos(Pg[K+1]);
            file_out.println(x+" "+y+" "+z);
	  } // end if
        } // end loop K
        file_out.println(" ");
      } // end loop J
    } // end loop I

    for(int J=0; J<Nt; J++)
    {
      for(int K=0; K<Np; K++)
      {
        for(int I=0; I<Nr-1; I++)
	{
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          file_out.println(x+" "+y+" "+z);
	  if(I+1<Nr)
	  {
	    x = Rg[I+1]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	    y = Rg[I+1]*Math.sin(Tg[J])*Math.sin(Pg[K]);
	    z = Rg[I+1]*Math.cos(Pg[K]);
            file_out.println(x+" "+y+" "+z);
	  } // end if
	} // end loop I
        file_out.println(" ");
      } // end loop K
    } // end loop J

    for(int K=0; K<Np; K++)
    {
      for(int I=0; I<Nr; I++)
      {
        for(int J=0; J<Nt-1; J++)
	{
	  x = Rg[I]*Math.cos(Tg[J])*Math.sin(Pg[K]);
	  y = Rg[I]*Math.sin(Tg[J])*Math.sin(Pg[K]);
          z = Rg[I]*Math.cos(Pg[K]);
          file_out.println(x+" "+y+" "+z);
          tgv = Tg[J];
	  if(J+1<Nt)
	  {
	    while(Tg[J+1]>tgv+0.1)
	    {
	      tgv = tgv + 0.1;
	      x = Rg[I]*Math.cos(tgv)*Math.sin(Pg[K]);
	      y = Rg[I]*Math.sin(tgv)*Math.sin(Pg[K]);
	      z = Rg[I]*Math.cos(Pg[K]);
              file_out.println(x+" "+y+" "+z);
	    } // end while
	    x = Rg[I]*Math.cos(Tg[J+1])*Math.sin(Pg[K]);
	    y = Rg[I]*Math.sin(Tg[J+1])*Math.sin(Pg[K]);
	    z = Rg[I]*Math.cos(Pg[K]);
            file_out.println(x+" "+y+" "+z);
	  } // end if
	} // end J
      } // end I
    } // end K

    file_out.close();
    System.out.println("wrote spherical_deriv.dat for gnuplot plotting ");
    System.out.println(" ");
  } // end write_file

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