// filename: align.java
// M.Lampton UCBSSL 2004
// mlampton@SSL.berkeley.edu
//
// Download and save this file in any directory.
// Name the file:  align.java
// To compile:     javac align.java
// To run it:      java align


import java.awt.*;
import java.awt.event.*;
import java.awt.geom.*;
import javax.swing.*;
import javax.swing.event.*;
import java.io.*;
import java.util.*; 
import java.text.DecimalFormat;

public class align
{
    public static void main(String[] args)
    {
        CFrame frame = new CFrame();
        frame.setSize(800, 650); 
        frame.setLocation(50, 50); 
        frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); 
        frame.show(); 
    }
}


class CFrame extends JFrame
{

    int[] data = new int[6]; 
    static final int DEFAULT_MICRONS = 100; 

    private JPanel radio;       /// class common is a convenience for addRadioButton()
    private ButtonGroup group;  /// class common is a convenience for addRadioButton()

    static final int NFIELDS = 4; // coordinate with NSTARS in User

    static JTextField myFields[] = new JTextField[NFIELDS];  // ..for updateField()
    static JTextField summary;  /// class common is a convenience for updateSummary()
    

    public CFrame()  // constructor for content frame
    {
        setTitle("     SNAP TELESCOPE SECONDARY MIRROR HEXAPOD CONTROLLER");
        Container contentPane = getContentPane();
        contentPane.setLayout(null);      // turn off layout managers!
       

        ////////////////////// top zone ////////////////////////////////

        JLabel t1 = new JLabel("Hexapod");
        t1.setBounds(50, 100, 100, 15); 
        contentPane.add(t1); 
        JLabel t2 = new JLabel("strut"); 
        t2.setBounds(50, 115, 100, 15); 
        contentPane.add(t2); 
        JLabel t3 = new JLabel("settings");
        t3.setBounds(50, 130, 100, 15); 
        contentPane.add(t3); 
        JLabel t4 = new JLabel("microns");
        t4.setBounds(50, 145, 100, 15); 
        contentPane.add(t4); 

        for (int i=0; i<6; i++)
        {
            final JSlider slider               // final for inner class below
              = new JSlider(JSlider.VERTICAL, -80, 80, 0);
            slider.setMajorTickSpacing(40);
            slider.setMinorTickSpacing(20);
            slider.setPaintTicks(true);
            slider.setPaintLabels(true); 
            slider.setSize(60,120); 
            slider.setBounds(120+80*i, 40, 50, 200); 
            contentPane.add(slider); 

            final JLabel myLabel               // final for inner class below
               = new JLabel(""+slider.getValue()); 
            data[i] = slider.getValue(); 
            myLabel.setBounds(130+80*i,20, 50, 20); 
            contentPane.add(myLabel); 

            final int ifinal = i;              // final for inner class below
            slider.addChangeListener(new ChangeListener(  ) 
            {
	        public void stateChanged(ChangeEvent e) 
                {
		   myLabel.setText(""+slider.getValue(  ));
                   data[ifinal] = slider.getValue(); 
	        }
            });
        }

        JButton myButton = new JButton("UPLINK"); 
        myButton.setActionCommand("EXECUTE"); 
        myButton.setSize(myButton.getPreferredSize()); 
        myButton.setBackground(Color.BLUE); 
        myButton.setForeground(Color.WHITE); 
        myButton.setBounds(650, 100, 90, 30); 
        contentPane.add(myButton); 

        JPanel buttonPanel = new JPanel(); 
        buttonPanel.setBounds(640, 90, 110, 50); 
        buttonPanel.setBackground(Color.YELLOW); 
        contentPane.add(buttonPanel); 

        //////////////// mid zone: announcement display //////////

        final JTextArea myText               // final for inner class below
           = new JTextArea(4, 36);           // 5 rows of 36 columns
        myText.setEditable(false);        
        JScrollPane sp = new JScrollPane(myText);  
        sp.setHorizontalScrollBarPolicy(ScrollPaneConstants.HORIZONTAL_SCROLLBAR_NEVER);    
        sp.setBounds(30, 270, 400, 90); 
        contentPane.add(sp); 


        /////////////// graphics displays ///////////

        radio = new JPanel(); 
        JLabel choose = new JLabel("Choose span of display fields:   "); 
        radio.add(choose); 
        group = new ButtonGroup();
        addRadioButton("500um", 500); 
        addRadioButton("200um", 200); 
        addRadioButton("100um", 100); 
        addRadioButton("50um",  50); 
        radio.setBounds(50, 365, 600, 25); 
        contentPane.add(radio); 


        PPanel[] mypanels = new PPanel[NFIELDS];    // names only
        for (int i=0; i<NFIELDS; i++)
        {
            mypanels[i] = new PPanel(i); 
            mypanels[i].setBounds(20+200*i, 395, 150, 150); 
            mypanels[i].setBackground(Color.WHITE); 
            contentPane.add(mypanels[i]); 
        }

        for (int i=0; i<NFIELDS; i++)
        {
            myFields[i] = new JTextField("no data", 20); 
            myFields[i].setBounds(20+200*i, 550, 150, 20); 
            contentPane.add(myFields[i]); 
        }

        summary = new JTextField("no summary data yet", 40); 
        summary.setBounds(300, 595, 200, 20); 
        contentPane.add(summary); 

        myButton.addActionListener( new ActionListener()
        {
            public void actionPerformed(ActionEvent ae)
            {
               myText.append("Uplinked:"+ni(6, data)+"\n"); 
               myText.append("Downlink received\n"); 
               User.recompute(data); 
               myText.append("Data have been processed\n"); 
               repaint(); 
            }
        });
    };

    public void addRadioButton(String name, int microns)
    {
        boolean selected = (microns == DEFAULT_MICRONS);
        JRadioButton button = new JRadioButton(name, selected); 
        group.add(button); 
        radio.add(button); 

        final int fmicrons = microns; 

        ActionListener listener = new ActionListener()  // H&C p.384
        {
           public void actionPerformed(ActionEvent ae)
           {
               User.setBoxSize(fmicrons); 
               repaint(); 
           }
        };
        button.addActionListener(listener); 
    }

    
    static void updateField(int whichfield, String s)
    // this is to allow User to update any field
    // to show results of calculations.
    {
        myFields[whichfield].setText(s); 
    }

    static void updateSummary(String s)
    // allows diagnostics to be posted on screen
    {
       summary.setText(s); 
    }


    static String ni(int n, int[] j)
    {
        String s = ""; 
        for (int i=0; i<n; i++)
          s += fwi(j[i], 6); 
        return s;
    }

    static String fwi(int n, int w)
    // converts an int to a string with given width.
    {
        String s = Integer.toString(n); 
        while (s.length() < w)
           s = " " + s; 
        return s;
    }
}







class PPanel extends JPanel  // integer plot panel
{
    private int mynumber; 

    public PPanel(int i)     // constructor! never "void"
    {
        mynumber = i; 
    }
    
    public void paintComponent(Graphics g)
    {
        super.paintComponent(g); 
        setBackground(Color.WHITE);      // set "this" panel background
        User.anyDraw(g, mynumber); 
    }

    ///////////// plotting symbols ///////////////

    static void putPixel(Graphics g, int i, int j)
    {
         g.drawLine(i, j, i, j); 
    }

    static void put2x2(Graphics g, int i, int j)
    {
         g.drawLine(i, j, i+1, j); 
         g.drawLine(i, j+1, i+1, j+1); 
    }

    static void putNxN(Graphics g, int i, int j, int s)
    // plots a solid square, center = (i,j)
    // s = size of symbol in pixels = 1, 2, ...
    {
        int start = -(s-1)/2; 
        int stop = s/2; 
        for (int k=start; k<=stop; k++)
          g.drawLine(i+start, j+k, i+stop, j+k); 
    }

    static void putPlus(Graphics g, int i, int j, int r)
    /// r = radius of symbol, in pixels
    {
        g.drawLine(i, j-r, i, j+r); 
        g.drawLine(i-r, j, i+r, j); 
    }

    static void putEx(Graphics g, int i, int j, int r)
    /// r = radius of symbol, in pixels
    {
        g.drawLine(i-r, j-r, i+r, j+r); 
        g.drawLine(i-r, j+r, i+r, j-r); 
    }
}





 
class User // performs user artwork using integer PPanel methods
{
    static final int NSTARS = 4; 
    static final double[][] uv = 
       {{0.01, 0.0}, {0, 0.012}, {-0.009, 0},{0, -0.008}}; 

    static double[][][] xylist = new double[NSTARS][331][2]; 
    static int[] xycount = new int[NSTARS]; 
    static double[][] centers = new double[NSTARS][2]; 

    static int iLeft = 0;
    static int iRight = 150; 
    static int jTop = 0; 
    static int jBot = 150; 

    static double boxsize = 1E-4;  // must coordinate with DEFAULT_BOXSIZE

    static double xLeft = -1; 
    static double xRight = 1; 
    static double yBot = -1; 
    static double yTop = 1;  

    static public void setBoxSize(double microns)
    {
        boxsize = 1E-6 * microns; 
    }

    static int geti(double x)
    {
       return (int)(iLeft + (x-xLeft)*(iRight-iLeft)/(xRight-xLeft)); 
    }

    static int getj(double y)
    {
       return (int)(jBot + (y-yBot)*(jTop-jBot)/(yTop-yBot)); 
    }

    static double sqr(double z)
    {
       return z*z; 
    }


    static void recompute(int[] data)
    {
         RT.setSurfs(); 

         double[] vec = new double[6]; 
         vec[0] = 1E-6*data[0]; 
         vec[1] = 1E-6*data[1]; 
         vec[2] = 1E-6*data[2]; 
         vec[3] = 1E-6*data[3]; 
         vec[4] = 1E-6*data[4]; 
         vec[5] = 1E-6*data[5];

         Stewart.convert(vec);          

         RT.surfs[2][RT.OX]     = vec[0];
         RT.surfs[2][RT.OY]     = vec[1];
         RT.surfs[2][RT.OZ]     = vec[2] - 2.000;
         RT.surfs[2][RT.OTILT]  = 10*vec[3]; 
         RT.surfs[2][RT.OPITCH] = 10*vec[4];
         RT.surfs[2][RT.OROLL]  = 10*vec[5];
 
         RT.setEulers(); 

         double[] rms = new double[NSTARS];  /// coordinate with NFIELDS

         for (int m=0; m<NSTARS; m++)
         {
             rms[m] = 0.0; 
             xycount[m] = 0; 
             for (int k=0; k<331; k++)
               if (RT.getxy(k, uv[m][0], uv[m][1], xylist[m][xycount[m]]))
                 xycount[m]++; 

             System.out.println("xycount = "+RT.fwi(xycount[m], 6)); 

             centers[m][0] = centers[m][1] = 0.0; 

             if (xycount[m] > 2)
             {
                 for (int k=0; k<xycount[m]; k++)
                 {
                     centers[m][0] += xylist[m][k][0];
                     centers[m][1] += xylist[m][k][1]; 
                 }
                 centers[m][0] /= xycount[m]; 
                 centers[m][1] /= xycount[m];

                 System.out.println("Centers =" + RT.nd(2, centers[m])); 

                 double x2sum=0, y2sum=0; 
                 for (int k=0; k<xycount[m]; k++)
                 {
                     x2sum += sqr(xylist[m][k][0]); 
                     y2sum += sqr(xylist[m][k][1]); 
                 }
                 double xrms = Math.sqrt((x2sum-xycount[m]*sqr(centers[m][0]))/xycount[m]); 
                 double yrms = Math.sqrt((y2sum-xycount[m]*sqr(centers[m][1]))/xycount[m]);
                 rms[m] = Math.sqrt(0.5*(sqr(xrms) + sqr(yrms))); 

                 // post this on the display

                 String s = "rms microns ="+ RT.fwd(1E6*rms[m], 12, 3); 
                 CFrame.updateField(m, s);   
            }
        } 

        ///////// post the overall summary /////////////////////
        double overall = 0.0; 
        for (int m=0; m<NSTARS; m++)
          overall += sqr(rms[m]); 

        overall = 1E6*Math.sqrt(overall/NSTARS); 
        String t = "Overall RMS microns = " + RT.fwd(overall, 12, 3); 
        CFrame.updateSummary(t); 
    }

    static void anyDraw(Graphics g, int i)
    {
        xLeft  = centers[i][0] - 0.5*boxsize; 
        xRight = centers[i][0] + 0.5*boxsize;
        yTop   = centers[i][1] + 0.5*boxsize;
        yBot   = centers[i][1] - 0.5*boxsize;  
 
        for (int k=0; k<xycount[i]; k++)
          PPanel.put2x2(g, geti(xylist[i][k][0]), getj(xylist[i][k][1])); 
    }
}






class Stewart
{
    static final double[][] linkstolab = // [row][col]
    {{ +0.66667, -0.66667, -0.33333, +0.33333, -0.33333, +0.33333},
     {  0.00000,  0.00000, -0.57737, +0.57737, +0.57737, -0.57737},
     { +0.19246, +0.19246, +0.19246, +0.19246, +0.19246, +0.19246},
     { +0.38491, +0.38491, -0.19246, -0.19246, -0.19246, -0.19246},
     {  0.00000,  0.00000, -0.33333, -0.33333, +0.33333, +0.33333},
     { -0.33333, +0.33333, -0.33333, +0.33333, -0.33333, +0.33333}}; 

    static void convert(double[] vec)
    {

         double[] result = new double[6]; 
         for (int i=0; i<6; i++)
         {
            result[i] = 0.0; 
            for (int k=0; k<6; k++)
              result[i] += linkstolab[i][k] * vec[k]; 
         }
         for (int i=0; i<6; i++)
           vec[i] = result[i]; 
    }
}
           
            




/********* class RT: ray tracing ************

        origin = vertex of primary mirror
        X, Y perpendicular to axis
        Z = axis
        all dimensions are in meters,
 
 M.Lampton STELLAR SOFTWARE  (c) 2003:  info@stellarsoftware.com
 Class RT is all static methods, no instantiation.
 Be sure to setSurfs() before tracing!
 This class includes a ray generator, raynumber 0...NRAYS
 The trace command is getxy(int raynumber, double xyresult[2])

*/

class RT
{
    static final boolean DEBUG = false; 

    ////////// RAY TRACING CODE STARTS HERE ///////////////

    static final double TINY  =1E-99;
    static final double TOL   =5e-14;
    static final double INF   =9e+99;
    static final double DELTA =1E-6; 
    static final double PI    =3.14159265358979324;
    static final double TWOPI =6.28318530717958648;
    static final double DEGRAD=57.2957795130823214;
    static final double LN10  =2.30258509299404568;
    static final double ROOT2 =1.41421356237309505;

    ///////// OPTICS TABLE INDEXING /////////////////////////

    static final int OX=0;             // LabX vertex
    static final int OY=1;             // LabY vertex
    static final int OZ=2;             // LabZ vertex
    static final int OTILT=3;          // tilt around Lab X, degrees
    static final int OPITCH=4;         // pitch around tilted Y, degrees
    static final int OROLL=5;          // roll around tilted & pitched Z
    static final int OTYPE=6;          // posvalues=refrIndex, neg=types
    static final int OFORM=7;          // 0.0=elliptical, 1.0=rectangle
    static final int OSTATUS=8;        // unused
    static final int OCURVE=12;        // 1/VertexRadius
    static final int OASPHER=13;       // 0=sphere, -1=parabola
    static final int OGX=14;           // grating grooveFreq,x
    static final int OGY=15;           // grating grooveFreq,y
    static final int OORDER=16;        // diffraction order: +-1, +-2,..
    static final int OSPARE=17;        // spare
    static final int OODIAX=18;        // outer diam, X
    static final int OODIAY=19;        // outer diam, Y
    static final int OIDIAX=20;        // inner diam, X
    static final int OIDIAY=21;        // inner diam, Y
    static final int OVLS1=22;
    static final int OVLS2=23;
    static final int OVLS3=24;
    static final int OVLS4=25;
    static final int OHOEX1=26;
    static final int OHOEY1=27;
    static final int OHOEZ1=28;
    static final int OHOEX2=29;
    static final int OHOEY2=30;
    static final int OHOEZ2=31;
    static final int OHOELAM=32;       // end of user parms
    static final int OE11=33;          // rotation matrix elements
    static final int OE12=34; 
    static final int OE13=35; 
    static final int OE21=36; 
    static final int OE22=37; 
    static final int OE23=38;
    static final int OE31=39;
    static final int OE32=40; 
    static final int OE33=41; 
    static final int ONPARMS=42; 

    //////// optics actions code numbers //////////////

    static final double OTMIRROR=-12.0;
    static final double OTIRIS=-13.0;

    // for reflection grating, just put grooves on a mirror 
    // for transmission grating, just put grooves on a lens

    ////////// optic element boundary shape codes //////

    static final double OFELLIP=0.0; 
    static final double OFRECT=1.0; 

    ////////// failure code numbers ///////////////

    static final int OK=0; 
    static final int PROPMISS=1; 
    static final int PROPBACK=2; 
    static final int PROPFAIL=3; 
    static final int IRISINNER=4; 
    static final int IRISOUTER=5; 
    static final int REDINRFAIL=6; 
    static final int REDGRFAIL=7; 
    static final int REDORDFAIL=8;
    static final int OFAIL=999; 

    /////////////// RAY VECTOR INDEXING ////////////////////

    static final int RX=0;
    static final int RY=1;
    static final int RZ=2; 
    static final int RU=3; 
    static final int RV=4; 
    static final int RW=5; 
    static final int RLAM=6; 
    static final int RNPARMS=7; 


    //////////// USER METHODS START HERE /////////////////////

    static public boolean getxy(int whichray, double u, double v, double[] xy)
    // Given whichray=0..330 of the iris group, Given U and V, 
    // runs one ray and returns its final focal plane x,y
    {
        double[] ray = new double[RNPARMS]; 
        ray[RU] = u;
        ray[RV] = v; 
        posWnormalize(ray);  // set ray[RW]
        build331Ray(ray, whichray); 
        int failcode = runray(ray); 
        xy[0] = ray[RX]; 
        xy[1] = ray[RY]; 
        return (failcode == 0);  
    }

    public static void build331Ray(double[] ray, int k)
    // Given a ray with RU, RV, RW, this routine
    // supplies RX, RY, RZ for one ray
    // Given index k=0...330 controls hex location within pupil
    // Note: hexagonal numbers 1, 2, 3, ...10 are 7, 19, 37, ...331.
    {
        double RADIUS = 1.0;              // radius of pupil
        int ORDER = 10;                   // 10 circles, 331 rays
        k = k % 331;                      // safety valve
        ray[RX] = 0.0;
        ray[RY] = 0.0; 
        ray[RZ] = -0.12; 

        int ktest = 0; 
        for (int j=ORDER; j>0; j--)   
        {
            for (int i=0; i<6*j; i++)    // do circle "j", or nothing if j==0.
            {
                if (k==ktest)
                {
                    double theta = 360.0/(6*j);
                    double thetaZero = isodd(j+ORDER) ? 0.5*theta : 0.0; 
                    ray[RX] = RADIUS*j*cosd(thetaZero + i*theta)/ORDER;
                    ray[RY] = RADIUS*j*sind(thetaZero + i*theta)/ORDER;
                }
                ktest++; 
            }
        }
        ray[RLAM] = 0.600; 

        if (k==0)
          System.out.println(RT.nd(6, ray)); 
        return; 
    }


    ////////////// optical surface definitions //////////////

    static final int ONSURFS = 5;     // surface zero is source 
    static double[][] surfs  = new double[ONSURFS+1][ONPARMS]; 

    static void setSurfs() //// global surfs 0, 1, 2, ...
    // call this once, before tracing, to set up surfs[][].
    // This is TMA63, as per SPIE 4849 p.215 Table 1.
    {
       for (int j=0; j<=ONSURFS; j++)
         for (int i=0; i<ONPARMS; i++)
            surfs[j][i] = 0.0; 

       int pri = 1; 
       surfs[pri][OZ] = 0.0; 
       surfs[pri][OCURVE] = -0.2037466;
       surfs[pri][OASPHER]= -0.981128; 
       surfs[pri][OTYPE] = OTMIRROR; 

       int sec = 2; 
       surfs[sec][OZ] = -2.0;
       surfs[sec][OCURVE] = -0.9099607;
       surfs[sec][OASPHER] = -1.847493; 
       surfs[sec][OTYPE] = OTMIRROR; 

       int flat = 3; 
       surfs[flat][OZ] = +0.91; 
       surfs[flat][OPITCH] = 45.0; 
       surfs[flat][OTYPE] = OTMIRROR; 

       int tert = 4; 
       surfs[tert][OZ] = 0.91; 
       surfs[tert][OX] = -0.87; 
       surfs[tert][OPITCH] = -90.0; 
       surfs[tert][OCURVE] = -0.7112388; 
       surfs[tert][OASPHER] = -0.599000; 
       surfs[tert][OTYPE] = OTMIRROR; 

       int fp = 5; 
       surfs[fp][OZ] = +0.91; 
       surfs[fp][OX] = +0.9; 
       surfs[fp][OPITCH] = 90.0; 

       setEulers(); 
    }



    ///////////////////////// ray tracer runray() ///////////////////////

         
    static int runray(double ray[])
    // Runs one ray through an optical system: global surfs[][]
    // returns a failurecode = surface where it stopped.
    // if failcode = 0 then no failure, went thru all.
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    {
        int errcode, tgt; 
        if (DEBUG) System.out.println("\n  Starting Ray"); 
        for (tgt=1; tgt<=ONSURFS; tgt++)  // tgt=1 is first target surface
        {
           if (DEBUG) System.out.println("Runray target="+fwi(tgt,3)); 
           if (DEBUG) System.out.println("LAB"+nd(6,ray)); 

           labtovx(ray, ray, surfs[tgt]);  
           if (DEBUG) System.out.println("VX="+nd(6,ray)); 

           errcode = intercept(ray, surfs[tgt]); 
           if (errcode != OK)
           {
              System.out.println("ERROR: intercept code=,  tgt=" + twoint(errcode,tgt)); 
              return tgt; 
           }
           if (DEBUG) System.out.println("VX="+nd(6,ray)); 

           if (tgt<ONSURFS)  // redirect the ray
           {
               double tnext = surfs[tgt+1][OTYPE];
               errcode = redirect(ray, surfs[tgt], tnext); 
               if (errcode != OK)
               {
                  System.out.println("ERROR: redirect code,tgt=" + twoint(errcode,tgt)); 
                  return tgt; 
               }
               if (DEBUG) System.out.println("VX="+nd(6,ray)); 
               vxtolab(ray, ray, surfs[tgt]); 
               if (DEBUG) System.out.println("LAB"+nd(6,ray)); 
           }
           //////// final surface is left in local detector frame
        }
        return 0;   
    }



    ///////////////////// coordinate utilities ////////////////////

    static void labtovx(double rLab[], double rVertex[], double surf[])
    // converts a lab-origin ray to a vertex-origin ray, one surface
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003
    {
        double a, b, c;
        a = rLab[RX] - surf[OX];
        b = rLab[RY] - surf[OY];
        c = rLab[RZ] - surf[OZ];
        rVertex[RX] = surf[OE11]*a + surf[OE12]*b + surf[OE13]*c;
        rVertex[RY] = surf[OE21]*a + surf[OE22]*b + surf[OE23]*c;
        rVertex[RZ] = surf[OE31]*a + surf[OE32]*b + surf[OE33]*c;

        a = rLab[RU];
        b = rLab[RV];
        c = rLab[RW];
        rVertex[RU] = surf[OE11]*a + surf[OE12]*b + surf[OE13]*c;
        rVertex[RV] = surf[OE21]*a + surf[OE22]*b + surf[OE23]*c;
        rVertex[RW] = surf[OE31]*a + surf[OE32]*b + surf[OE33]*c;
    }


    static void vxtolab(double rVertex[], double rLab[], double surf[])
    // converts a vertex-origin ray to a lab-origin ray, one surface
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003
    {
        double a = rVertex[RU];
        double b = rVertex[RV];
        double c = rVertex[RW];
        rLab[RU] = surf[OE11]*a + surf[OE21]*b + surf[OE31]*c;
        rLab[RV] = surf[OE12]*a + surf[OE22]*b + surf[OE32]*c;
        rLab[RW] = surf[OE13]*a + surf[OE23]*b + surf[OE33]*c;
        a = rVertex[RX];
        b = rVertex[RY];
        c = rVertex[RZ];
        rLab[RX] = surf[OE11]*a + surf[OE21]*b + surf[OE31]*c;
        rLab[RY] = surf[OE12]*a + surf[OE22]*b + surf[OE32]*c;
        rLab[RZ] = surf[OE13]*a + surf[OE23]*b + surf[OE33]*c;
        rLab[RX] += surf[OX];
        rLab[RY] += surf[OY];
        rLab[RZ] += surf[OZ];
    }

    static void setEulers()
    // sets up Euler matrix coefs for surfs[][].
    // sequence is tilt(x), pitch(y), roll(z).
    // call this after setting up surface, before ray tracing. 
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    {
       double ct, st, cp, sp, cr, sr; 
       for (int j=1; j<=ONSURFS; j++)
       {
           ct = cosd(surfs[j][OTILT]);           // degrees
           st = sind(surfs[j][OTILT]);           // degrees
           cp = cosd(surfs[j][OPITCH]);          // degrees
           sp = sind(surfs[j][OPITCH]);          // degrees
           cr = cosd(surfs[j][OROLL]);           // degrees
           sr = sind(surfs[j][OROLL]);           // degrees
           surfs[j][OE11] = cr*cp;               // x <-> X
           surfs[j][OE12] = cr*sp*st + sr*ct;    // x <-> Y
           surfs[j][OE13] = -cr*sp*ct + sr*st;   // x <-> Z
           surfs[j][OE21] = -sr*cp;              // y <-> X
           surfs[j][OE22] = cr*ct - sr*sp*st;    // y <-> Y
           surfs[j][OE23] = sr*sp*ct + cr*st;    // y <-> Z
           surfs[j][OE31] = sp;                  // z <-> X
           surfs[j][OE32] = -cp*st;              // z <-> Y
           surfs[j][OE33] = cp*ct;               // z <-> Z
       }
    }    

    ///////////// optical utilities /////////////////////////////

    static int intercept(double ray[], double surf[])
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    {
        if (sqr(surf[OCURVE]) < TOL)
          return planesolve(ray, surf);
        return quadsolve(ray, surf); 
    }


    static int redirect(double ray[], double surf[], double tnext)
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    {
        if (surf[OTYPE]==OTIRIS)      // iris
          return iris(ray, surf);

        boolean groovy = sqr(surf[OGX]) + sqr(surf[OGY]) > TINY; 

        if (surf[OTYPE]==OTMIRROR)    // mirror
          if (groovy) 
            return rgrating(ray, surf);
          else
            return mirror(ray, surf);

        if (groovy)                   // lens
          return tgrating(ray, surf, tnext);
        else
          return snell(ray, surf, tnext); 
    }


    static int rgrating(double ray[], double surf[])
    {
        return OFAIL; 
    }

    static int tgrating(double ray[], double surf[], double t)
    {
        return OFAIL; 
    }

    static int snell(double ray[], double surf[], double t)
    {
        return OFAIL; 
    }

    static int planesolve(double ray[], double surf[])
    // if valid plane intercept, ray[] gets updated, return OK
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    // NEEDED: diameter check
    {
        double d;

        if (Math.abs(ray[RZ]) < TOL)
          return OK; 

        if (Math.abs(ray[RW]) < TOL)
          return PROPMISS;

        d = -ray[RZ] / ray[RW];
        if (d < 0.0)                       
          return PROPBACK;                   

        ray[RX] += d * ray[RU]; 
        ray[RY] += d * ray[RV]; 
        ray[RZ] += d * ray[RW]; 
        return OK;
    }


    static int quadsolve(double ray[], double surf[])
    // if valid quadratic surface intercept, ray[] updated, return OK
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    // NEEDED: diameter testing
    {
        double a, b, c, k, s, g;
        double rx, ry, rz, ru, rv, rw;

        if (sqr(surf[OCURVE]) < TOL)
          System.out.println("ERROR: Quadsolve was called for a plane intercept!");

        // shorthand variables...
        k = surf[OCURVE];         // curvature
        s = surf[OASPHER]+1.0;    // shape
        g = 1.0;                  // X-stretch is unused
        rx = ray[RX];
        ry = ray[RY];
        rz = ray[RZ];
        ru = ray[RU];
        rv = ray[RV];
        rw = ray[RW];
        a = k * (sqr(g*ru) + sqr(rv) + s*sqr(rw));
        b = 2.0 * k * (g*g*rx*ru + ry*rv + s*rz*rw) - 2.0*rw;
        c = k * (sqr(g*rx) + sqr(ry) + s*sqr(rz)) - 2.0*rz;

        double[] d = new double[2]; 
        int nposroots = getposroots(a, b, c, d); 
        if (DEBUG) System.out.println("getposroots() returns nroots = "+fwi(nposroots,3)); 
        switch (nposroots)
        {
            case 0: 
              return PROPFAIL;

            case 1:  // here I need to test iris() before agreeing
              ray[RX] += d[0] * ray[RU]; 
              ray[RY] += d[0] * ray[RV]; 
              ray[RZ] += d[0] * ray[RW]; 
              return OK;

            case 2:  // choose root nearer vertex
              double azvx0 = Math.abs(ray[RZ] + d[0]*ray[RW]);      
              double azvx1 = Math.abs(ray[RZ] + d[1]*ray[RW]);
              double dd = (azvx0 < azvx1) ? d[0] : d[1]; 
              ray[RX] += dd * ray[RU]; 
              ray[RY] += dd * ray[RV]; 
              ray[RZ] += dd * ray[RW]; 
              return OK; 
        }
        return PROPFAIL; 
    }


    static int iris(double ray[], double surf[])
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    /// ..a better way is to avoid ratios, and always compare differences:
    /// diff[i] = ray - iris;
    {
        boolean square = (surf[OFORM] == OFRECT);
 
        // do the inner edge
        if ((surf[OIDIAX] > TINY) && (surf[OIDIAY] > TINY))
        {
          if (square &&
             (max(sqr((ray[RX])/surf[OIDIAX]), sqr((ray[RY])/surf[OIDIAY])) < 0.25))
          return IRISINNER;
          if ((!square) &&
             (sqr((ray[RX])/surf[OIDIAX]) + sqr((ray[RY])/surf[OIDIAY]) < 0.25))
          return IRISINNER;
        }

        // do the outer edge
        if ((surf[OODIAX] > TINY) && (surf[OODIAY] > TINY))
        {
          if (square &&
             (max(sqr((ray[RX])/surf[OODIAX]), sqr((ray[RY])/surf[OODIAY])) > 0.25))
          return IRISOUTER;

          if (!square &&
           (sqr((ray[RX])/surf[OODIAX]) + sqr((ray[RY])/surf[OODIAY]) > 0.25))
          return IRISOUTER;
        }
        return OK; 
    }





    static int mirror(double ray[], double surf[])
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    // method: r = i - 2 (i dot n) n
    // note this is quadratic in n, hence independent of its sign
    {
        double[] norm = new double[6];
        getnormal(ray, surf, norm); 
        double dotin = ray[RU]*norm[RU] + ray[RV]*norm[RV] + ray[RW]*norm[RW];
        ray[RU] -= 2.0 * dotin * norm[RU];
        ray[RV] -= 2.0 * dotin * norm[RV];
        ray[RW] -= 2.0 * dotin * norm[RW];
        normalize(ray);   // should be unnecessary!
        return OK;
    }




    static void getnormal(double ray[], double surf[], double normal[])
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    {
        /// analytic method: Spencer & Murty eq. 18, conics only
        /// first, verify that func() returns approx zero:
        double derr = getfunc(ray, surf); 
        if (DEBUG) System.out.println("Getnormal sees fn_val="+fwe(derr)); 

        // now get analytic gradient
        double r2 = sqr(ray[RX]) + sqr(ray[RY]); 
        double arg = 1.0 - (surf[OASPHER]+1.0)*sqr(surf[OCURVE])*r2; 
        if (arg > 0.0)   // valid surface...
        {
           double coef = surf[OCURVE]/Math.sqrt(arg); 
           normal[RU] = -ray[RX]*coef; 
           normal[RV] = -ray[RY]*coef; 
           normal[RW] = 1.0; 
           normalize(normal); 
        }
        else            // return flange normal...
        {
            normal[RU] = 0.0; 
            normal[RV] = 0.0; 
            normal[RW] = 1.0; 
        };

        if (DEBUG)  // compare with finite difference method
        {
            double delta = DELTA / max(1.0, Math.abs(surf[OCURVE])); 
            double[] fdnorm = new double[6]; 
            double[] rayp = new double[6]; 
            double[] raym = new double[6]; 
            for (int i=0; i<6; i++)
              rayp[i] = raym[i] = ray[i]; 
            rayp[RX] += delta; 
            raym[RX] -= delta; 
            double fxp = getfunc(rayp, surf); 
            double fxm = getfunc(raym, surf); 
            rayp[RX] = raym[RX] = ray[RX]; 
            rayp[RY] += delta; 
            raym[RY] -= delta; 
            double fyp = getfunc(rayp, surf); 
            double fym = getfunc(raym, surf); 
            rayp[RY] = raym[RY] = ray[RY]; 
            rayp[RZ] += delta; 
            raym[RZ] -= delta; 
            double fzp = getfunc(rayp, surf); 
            double fzm = getfunc(raym, surf); 
            fdnorm[RU] = fxp-fxm;
            fdnorm[RV] = fyp-fym;
            fdnorm[RW] = fzp-fzm; 
            normalize(fdnorm); 
        
            double fderr = Math.sqrt(sqr(normal[RU]-fdnorm[RU])
                                   + sqr(normal[RV]-fdnorm[RV])
                                   + sqr(normal[RW]-fdnorm[RW])); 
            System.out.println("Getnormal gets fd_err="+fwe(fderr));
        }
    }


    static double getfunc(double ray[], double surf[])
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    // returns z-zsurf(x,y) for flanged quadratic.
    // ray[] lies on surface if getfunc()==0. 
    {
        double r2 = sqr(ray[RX]) + sqr(ray[RY]); 
        double num = surf[OCURVE] * r2; 
        double arg = 1.0 - (surf[OASPHER]+1.0)*sqr(surf[OCURVE])*r2;
        arg = max(arg, 0.0); 
        double denom = 1.0 + Math.sqrt(arg); 
        double zsurf = num/denom; 
        return ray[RZ] - zsurf;
    }


    static void normalize(double norm[])
    {
        double r2 = sqr(norm[RU]) + sqr(norm[RV]) + sqr(norm[RW]); 
        if (r2 > TINY)
        {
            double r = Math.sqrt(r2); 
            norm[RU] /= r; 
            norm[RV] /= r; 
            norm[RW] /= r; 
            return; 
        }
        norm[RU] = 0.0; 
        norm[RV] = 0.0; 
        norm[RW] = 1.0; 
        System.out.println("ERROR: Normalize has substituted default unit vector"); 
    }


    static void posWnormalize(double ray[])
    {
        double w2 = 1.0 - sqr(ray[RU]) - sqr(ray[RV]); 
        if (w2 >= 0.0)
        {
           ray[RW] = Math.sqrt(w2); 
           return; 
        }
        normalize(ray); 
        System.out.println("ERROR: posWnormalize has had to call normalize"); 
    }
           

    static void negWnormalize(double[] ray)
    {
        double w2 = 1.0 - sqr(ray[RU]) - sqr(ray[RV]); 
        if (w2 >= 0.0)
        {
           ray[RW] = -Math.sqrt(w2); 
           return; 
        }
        normalize(ray); 
        System.out.println("ERROR: negWnormalize has had to call normalize"); 
    }
           
    ///////////////////// math utilities ////////////////////

    static double sqr(double x) 
    {
        return x*x;
    }

    static double max(double a, double b)
    {
        if (a>b) return a;
        return b;
    }

    static double min(double a, double b)
    {
        if (a<b) return a;
        return b; 
    }


    static double cosd(double x)
    {
        return Math.cos(x/DEGRAD); 
    }

    static double sind(double x)
    {
        return Math.sin(x/DEGRAD); 
    }

    static double tand(double x)
    {
        return Math.tan(x/DEGRAD); 
    }
    
    
    static String fwi(int n, int w)
    // converts an int to a string with given width.
    {
        String s = Integer.toString(n); 
        while (s.length() < w)
           s = " " + s; 
        return s;
    }


    static String twoint(int i, int j)
    // converts two small ints to a string of width 8
    {
       return fwi(i,4) + fwi(j,4); 
    }


    static String nd(int n, double[] v)
    // converts a vector of n doubles into a string
    {
       StringBuffer sb = new StringBuffer(100);
       for (int i=0; i<n; i++)
          sb.append(fwd(v[i],12,6)); 
       String s = sb.toString();
       return s; 
    }


    static String fwd(double x, int w, int d)
    // converts a double to a string with given width and decimals.
    {
        java.text.DecimalFormat df = new DecimalFormat();
        df.setMaximumFractionDigits(d); 
        df.setMinimumFractionDigits(d);
        df.setGroupingUsed(false);   
        String s = df.format(x); 
        while (s.length() < w)
          s = " " + s;  
        if (s.length() > w)
        {
            s = "";
            for (int i=0; i<w; i++)
              s = s + "-";
        }  
        return s; 
    }

    static String fwe(double x)
    // converts a double into an E notation string
    {
        java.text.NumberFormat ef = new DecimalFormat("0.00E0"); 
        String s = ef.format(x);    
        while (s.length() < 10)
          s = " " + s; 
        return s; 
    }

    static String fws(String s, int len)
    // enlarges a String so that it occupies len characters
    {
        while (s.length() < len)
          s = s + " ";
        return s; 
    }

    static double log10(double arg)
    {
        if (arg>0.0)
          return Math.log(arg)/LN10;
        return -999.9;
    }

    static double dex(double arg)
    {
        if (arg > 300.0)
           arg = 300.0; 
        if (arg < -300.0)
           arg = -300.0; 
        return Math.pow(10.0, arg); 
    }

    static boolean isodd(int i)
    {
        return (i % 2) == 1;
    }

    static void beep()
    {
        char c = 7; 
        String s = ""+c; 
        System.out.print(s); 
    }

    static int getposroots(double a, double b, double c, double d[])
    // Returns 0, 1,or 2 positive real roots of ax2+bx+c
    // Method derives from Press et al "Numerical Recipes" 5.5.5
    // M.Lampton STELLAR SOFTWARE (C) 1989, 2003 
    {
        double discrim, q, r1, r2;
        d[0] = d[1] = -999.0;          // errorcode

        if (Math.abs(a) < TOL)         // LINEAR FUNCTION
          if (Math.abs(b) > TOL)       // one root
            if (-c/b >= 0.0)           // nonnegative?
            {
               d[0] = d[1] = -c/b;
               if (DEBUG) System.out.println("Posroots exit 1,"+nd(1,d)); 
               return 1;               // one positive root
            }
            else
            {
               if (DEBUG) System.out.println("Posroots exit 2, one neg root"); 
               return 0;               // one nonpositive root
            }
          else 
          {
              if (DEBUG) System.out.println("Posroots exit 3, no roots"); 
              return 0;               // no roots
          };

        if (Math.abs(c) < TOL)         // LINEAR FUNCTION
          if (-b/a >= 0.0)             // nonnegative?
          {
              d[0] = d[1] = -b/a;
              if (DEBUG) System.out.println("Posroots exit 4,"+nd(1,d)); 
              return 1;                // one positive root
          }
          else 
          {
              if (DEBUG) System.out.println("Posroots exit 5, one neg root"); 
              return 0;                // one nonpositive root
          };

        discrim = b*b-4.0*a*c;         // QUADRATIC FUNCTION
        if (discrim < -TOL)            // no real roots
        {
            if (DEBUG) System.out.println("Posroots exit 6, no real roots"); 
            return 0;
        }

        if (discrim < TOL)             // degenerate pair at -b/2a
          if (-0.5*b/a > TOL)
          {
              d[0]=d[1]= -0.5*b/a;     // return positive pair as one
              if (DEBUG) System.out.println("Posroots exit 7,"+nd(1,d)); 
              return 1;
          }
          else 
          {
              if (DEBUG) System.out.println("Posroots exit 8, neg pair"); 
              return 0;                // nonpositive pair
          };

        if (b > 0.0)                   // two distinct roots
          q = -0.5*(b+Math.sqrt(discrim));
        else
          q = -0.5*(b-Math.sqrt(discrim));

        r1 = min(c/q, q/a);
        r2 = max(c/q, q/a);

        if (r1 >= 0.0)                 // is the smaller positive?
        {
            d[0] = r1;
            d[1] = r2;
            if (DEBUG) System.out.println("Posroots exit 9,"+nd(2,d)); 
            return 2;
        }
        if (r2 > TOL)                  // is the larger positive?
        {
            d[0] = d[1] = r2;
            if (DEBUG) System.out.println("Posroots exit 10,"+nd(1,d)); 
            return 1;
        }
        if (DEBUG) System.out.println("Posroots exit 11, two neg roots"); 
        return 0;                      // neither is positive
    }
}