/*

 ============================================================================
                   The Apache Software License, Version 1.1
 ============================================================================

 Copyright (C) 1999-2003 The Apache Software Foundation. All rights reserved.

 Redistribution and use in source and binary forms, with or without modifica-
 tion, are permitted provided that the following conditions are met:

 1. Redistributions of  source code must  retain the above copyright  notice,
    this list of conditions and the following disclaimer.

 2. Redistributions in binary form must reproduce the above copyright notice,
    this list of conditions and the following disclaimer in the documentation
    and/or other materials provided with the distribution.

 3. The end-user documentation included with the redistribution, if any, must
    include  the following  acknowledgment:  "This product includes  software
    developed  by the  Apache Software Foundation  (http://www.apache.org/)."
    Alternately, this  acknowledgment may  appear in the software itself,  if
    and wherever such third-party acknowledgments normally appear.

 4. The names "Batik" and  "Apache Software Foundation" must  not  be
    used to  endorse or promote  products derived from  this software without
    prior written permission. For written permission, please contact
    apache@apache.org.

 5. Products  derived from this software may not  be called "Apache", nor may
    "Apache" appear  in their name,  without prior written permission  of the
    Apache Software Foundation.

 THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED WARRANTIES,
 INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
 FITNESS  FOR A PARTICULAR  PURPOSE ARE  DISCLAIMED.  IN NO  EVENT SHALL  THE
 APACHE SOFTWARE  FOUNDATION  OR ITS CONTRIBUTORS  BE LIABLE FOR  ANY DIRECT,
 INDIRECT, INCIDENTAL, SPECIAL,  EXEMPLARY, OR CONSEQUENTIAL  DAMAGES (INCLU-
 DING, BUT NOT LIMITED TO, PROCUREMENT  OF SUBSTITUTE GOODS OR SERVICES; LOSS
 OF USE, DATA, OR  PROFITS; OR BUSINESS  INTERRUPTION)  HOWEVER CAUSED AND ON
 ANY  THEORY OF LIABILITY,  WHETHER  IN CONTRACT,  STRICT LIABILITY,  OR TORT
 (INCLUDING  NEGLIGENCE OR  OTHERWISE) ARISING IN  ANY WAY OUT OF THE  USE OF
 THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 This software  consists of voluntary contributions made  by many individuals
 on  behalf of the Apache Software  Foundation. For more  information on the
 Apache Software Foundation, please see <http://www.apache.org/>.

*/

package org.apache.batik.ext.awt;

import java.awt.Color;
import java.awt.Rectangle;
import java.awt.RenderingHints;
import java.awt.geom.AffineTransform;
import java.awt.geom.NoninvertibleTransformException;
import java.awt.geom.Rectangle2D;
import java.awt.image.ColorModel;

/**
 * Provides the actual implementation for the RadialGradientPaint.
 * This is where the pixel processing is done.  A RadialGradienPaint
 * only supports circular gradients, but it should be possible to scale
 * the circle to look approximately elliptical, by means of a
 * gradient transform passed into the RadialGradientPaint constructor.
 *
 * @author Nicholas Talian, Vincent Hardy, Jim Graham, Jerry Evans
 * @author <a href="mailto:vincent.hardy@eng.sun.com">Vincent Hardy</a>
 * @version $Id: RadialGradientPaintContext.java,v 1.10 2003/08/09 16:58:43 deweese Exp $
 *
 */
final class RadialGradientPaintContext extends MultipleGradientPaintContext {  
    
    /** True when (focus == center)  */
    private boolean isSimpleFocus = false;

    /** True when (cycleMethod == NO_CYCLE) */
    private boolean isNonCyclic = false;
       
    /** Radius of the outermost circle defining the 100% gradient stop. */
    private float radius;   
    
    /** Variables representing center and focus points. */
    private float centerX, centerY, focusX, focusY;     

    /** Radius of the gradient circle squared. */
    private float radiusSq; 
        
    /** Constant part of X, Y user space coordinates. */
    private float constA, constB;
       
    /** This value represents the solution when focusX == X.  It is called
     * trivial because it is easier to calculate than the general case.
     */
    private float trivial;       

    private static final int FIXED_POINT_IMPL = 1;
    private static final int DEFAULT_IMPL     = 2;
    private static final int ANTI_ALIAS_IMPL  = 3;

    private int fillMethod;
    
    /** Amount for offset when clamping focus. */
    private static final float SCALEBACK = .97f;
    
    /** 
     * Constructor for RadialGradientPaintContext.
     *
     *  @param cm {@link ColorModel} that receives
     *  the <code>Paint</code> data. This is used only as a hint.
     *
     *  @param deviceBounds the device space bounding box of the 
     *  graphics primitive being rendered
     *
     *  @param userBounds the user space bounding box of the 
     *  graphics primitive being rendered
     * 
     *  @param t the {@link AffineTransform} from user
     *  space into device space (gradientTransform should be 
     *  concatenated with this)    
     *
     *  @param hints the hints that the context object uses to choose
     *  between rendering alternatives
     *       
     *  @param cx the center point in user space of the circle defining 
     *  the gradient.  The last color of the gradient is mapped to the 
     *  perimeter of this circle X coordinate
     *
     *  @param cy the center point in user space of the circle defining 
     *  the gradient.  The last color of the gradient is mapped to the 
     *  perimeter of this circle Y coordinate
     *     
     *  @param r the radius of the circle defining the extents of the 
     *  color gradient
     *
     *  @param fx the point in user space to which the first color is mapped
     *  X coordinate
     *
     *  @param fy the point in user space to which the first color is mapped
     *  Y coordinate
     *
     *  @param fractions the fractions specifying the gradient distribution
     *
     *  @param colors the gradient colors
     *
     *  @param cycleMethod either NO_CYCLE, REFLECT, or REPEAT
     *
     *  @param colorSpace which colorspace to use for interpolation, 
     *  either SRGB or LINEAR_RGB
     *
     */
    public RadialGradientPaintContext(ColorModel cm,
                                      Rectangle deviceBounds,
                                      Rectangle2D userBounds,
                                      AffineTransform t,
                                      RenderingHints hints,
                                      float cx, float cy,
                                      float r,
                                      float fx, float fy,
                                      float[] fractions,
                                      Color[] colors,
                                      MultipleGradientPaint.CycleMethodEnum 
                                      cycleMethod,
                                      MultipleGradientPaint.ColorSpaceEnum 
                                      colorSpace)
        throws NoninvertibleTransformException
    {       	
        super(cm, deviceBounds, userBounds, t, hints, fractions, colors, 
              cycleMethod, colorSpace);

        //copy some parameters.
        centerX = cx;
        centerY = cy;	
        focusX = fx;
        focusY = fy;
        radius = r;

        this.isSimpleFocus = (focusX == centerX) && (focusY == centerY);
        this.isNonCyclic = (cycleMethod == RadialGradientPaint.NO_CYCLE);
	
        //for use in the quadractic equation
        radiusSq = radius * radius;

        float dX = focusX - centerX;
        float dY = focusY - centerY;

        double dist = Math.sqrt((dX * dX) + (dY * dY));

        //test if distance from focus to center is greater than the radius
        if (dist > radius* SCALEBACK) { //clamp focus to radius
          double angle = Math.atan2(dY, dX);
			
          //x = r cos theta, y = r sin theta
          focusX = (float)(SCALEBACK * radius * Math.cos(angle)) + centerX;
          
          focusY = (float)(SCALEBACK * radius * Math.sin(angle)) + centerY;
        }

        //calculate the solution to be used in the case where X == focusX
        //in cyclicCircularGradientFillRaster
        dX = focusX - centerX;
        trivial = (float)Math.sqrt(radiusSq - (dX * dX));

        // constant parts of X, Y user space coordinates 
        constA = a02 - centerX;
        constB = a12 - centerY;

        Object colorRend = hints.get(RenderingHints.KEY_COLOR_RENDERING);
        Object rend      = hints.get(RenderingHints.KEY_RENDERING);

        fillMethod = 0;

        if ((rend      == RenderingHints.VALUE_RENDER_QUALITY) ||
            (colorRend == RenderingHints.VALUE_COLOR_RENDER_QUALITY)) {
            // System.out.println("AAHints set: " + rend + ", " + colorRend);
            fillMethod = ANTI_ALIAS_IMPL;
        }

        if ((rend      == RenderingHints.VALUE_RENDER_SPEED) ||
            (colorRend == RenderingHints.VALUE_COLOR_RENDER_SPEED)) {
            // System.out.println("SPHints set: " + rend + ", " + colorRend);
            fillMethod = DEFAULT_IMPL;
        }

        // We are in the 'default' case, no hint or hint set to
        // DEFAULT values...
        if (fillMethod == 0) {
            // For now we will always use the 'default' impl if 
            // one is not specified.
            fillMethod = DEFAULT_IMPL;

            if (false) {
                // This could be used for a 'smart' choice in
                // the default case, if the gradient has obvious
                // discontinuites use AA, otherwise default
                if (hasDiscontinuity) {
                    fillMethod = ANTI_ALIAS_IMPL;
                } else {
                    fillMethod = DEFAULT_IMPL;
                }
            }
        }

        if ((fillMethod == DEFAULT_IMPL) &&
            (isSimpleFocus && isNonCyclic && isSimpleLookup)) {
            this.calculateFixedPointSqrtLookupTable();
            fillMethod = FIXED_POINT_IMPL;
        }
    }
    
    /**
     * Return a Raster containing the colors generated for the graphics
     * operation.
     * @param x,y,w,h The area in device space for which colors are
     * generated.
     */
    protected void fillRaster(int pixels[], int off, int adjust,
                              int x, int y, int w, int h) {
        switch(fillMethod) {
        case FIXED_POINT_IMPL:
            // System.out.println("Calling FP");
            fixedPointSimplestCaseNonCyclicFillRaster(pixels, off, adjust, x, 
                                                      y, w, h);
            break;
        case ANTI_ALIAS_IMPL:
            // System.out.println("Calling AA");
            antiAliasFillRaster(pixels, off, adjust, x, y, w, h);
            break;
        case DEFAULT_IMPL:
        default:
            // System.out.println("Calling Default");
            cyclicCircularGradientFillRaster(pixels, off, adjust, x, y, w, h);
        }
    }    
    
    /**
     * This code works in the simplest of cases, where the focus == center 
     * point, the gradient is noncyclic, and the gradient lookup method is 
     * fast (single array index, no conversion necessary).
     *
     */         
    private void fixedPointSimplestCaseNonCyclicFillRaster(int pixels[], 
                                                           int off,
                                                           int adjust, 
                                                           int x, int y, 
                                                           int w, int h) {
        float iSq=0;  // Square distance index
        final float indexFactor = fastGradientArraySize / radius;      

        //constant part of X and Y coordinates for the entire raster
        final float constX = (a00*x) + (a01*y) + constA;
        final float constY = (a10*x) + (a11*y) + constB;
        final float deltaX = indexFactor * a00; //incremental change in dX
        final float deltaY = indexFactor * a10; //incremental change in dY
        float dX, dY; //the current distance from center
        final int fixedArraySizeSq=
            (fastGradientArraySize * fastGradientArraySize);
        float g, gDelta, gDeltaDelta, temp; //gradient square value
        int gIndex; // integer number used to index gradient array
        int iSqInt; // Square distance index       		   
	
        int end, j; //indexing variables
        int indexer = off;//used to index pixels array

        temp        = ((deltaX * deltaX) + (deltaY * deltaY));
        gDeltaDelta = ((temp * 2));

        if (temp > fixedArraySizeSq) {
            // This combination of scale and circle radius means
            // essentially no pixels will be anything but the end
            // stop color.  This also avoids math problems.
            final int val = gradientOverflow;
            for(j = 0; j < h; j++){ //for every row
                //for every column (inner loop begins here)
                for (end = indexer+w; indexer < end; indexer++) 
                    pixels[indexer] = val;
                indexer += adjust;
            }
            return;
        }

        // For every point in the raster, calculate the color at that point
        for(j = 0; j < h; j++){ //for every row
            //x and y (in user space) of the first pixel of this row
            dX = indexFactor * ((a01*j) + constX);
            dY = indexFactor * ((a11*j) + constY);	   	   

            // these values below here allow for an incremental calculation
            // of dX^2 + dY^2 

            //initialize to be equal to distance squared
            g = (((dY * dY) + (dX * dX)) );
            gDelta =  (((((deltaY * dY) + (deltaX * dX))* 2) + 
                        temp));	 
	    
            //for every column (inner loop begins here)
            for (end = indexer+w; indexer < end; indexer++) {	       
                //determine the distance to the center
		
                //since this is a non cyclic fill raster, crop at "1" and 0
                if (g >= fixedArraySizeSq) {
                    pixels[indexer] = gradientOverflow;
                }
		
                // This should not happen as gIndex is a square
                // quantity. Code commented out on purpose, can't underflow.
                // else if (g < 0) {
                //    gIndex = 0;		    
                // }
		
                else {
                    iSq = (g * invSqStepFloat);
                    
                    iSqInt = (int)iSq; //chop off fractional part
                    iSq -= iSqInt;		    
                    gIndex = sqrtLutFixed[iSqInt];
                    gIndex += (int)(iSq * (sqrtLutFixed[iSqInt + 1]-gIndex));
                    pixels[indexer] = gradient[gIndex]; 
                }
		
				
                //incremental calculation
                g += gDelta;
                gDelta += gDeltaDelta;		
            }	  
            indexer += adjust;
        }
    }

    /** Length of a square distance intervale in the lookup table */
    private float invSqStepFloat; 
    
    /** Used to limit the size of the square root lookup table */
    private int MAX_PRECISION = 256;
    
    /** Square root lookup table */
    private int sqrtLutFixed[] = new int[MAX_PRECISION];
    
    /**
     * Build square root lookup table
     */       
    private void calculateFixedPointSqrtLookupTable() {	      
        float sqStepFloat;
        sqStepFloat = ((fastGradientArraySize  * fastGradientArraySize) 
                       / (MAX_PRECISION - 2));
	
        // The last two values are the same so that linear square root 
        // interpolation can happen on the maximum reachable element in the 
        // lookup table (precision-2)
        int i;
        for (i = 0; i < MAX_PRECISION - 1; i++) {
            sqrtLutFixed[i] = (int)(Math.sqrt(i*sqStepFloat));
        }
        sqrtLutFixed[i] = sqrtLutFixed[i-1];	
        invSqStepFloat = 1/sqStepFloat;
    }
    
    /** Fill the raster, cycling the gradient colors when a point falls outside
     *  of the perimeter of the 100% stop circle.          
     * 
     *  This calculation first computes the intersection point of the line
     *  from the focus through the current point in the raster, and the
     *  perimeter of the gradient circle.
     * 
     *  Then it determines the percentage distance of the current point along
     *  that line (focus is 0%, perimeter is 100%). 
     *
     *  Equation of a circle centered at (a,b) with radius r:
     *  (x-a)^2 + (y-b)^2 = r^2
     *  Equation of a line with slope m and y-intercept b
     *  y = mx + b
     *  replacing y in the cirlce equation and solving using the quadratic
     *  formula produces the following set of equations.  Constant factors have
     *  been extracted out of the inner loop.
     *
     */   
    private void cyclicCircularGradientFillRaster(int pixels[], int off, 
                                                  int adjust, 
                                                  int x, int y, 
                                                  int w, int h) {
        // Constant part of the C factor of the quadratic equation
        final double constC = 
            -(radiusSq) + (centerX * centerX) + (centerY * centerY);
        double A; //coefficient of the quadratic equation (Ax^2 + Bx + C = 0)
        double B; //coefficient of the quadratic equation
        double C; //coefficient of the quadratic equation
        double slope; //slope of the focus-perimeter line
        double yintcpt; //y-intercept of the focus-perimeter line
        double solutionX;//intersection with circle X coordinate
        double solutionY;//intersection with circle Y coordinate       
       	final float constX = (a00*x) + (a01*y) + a02;//const part of X coord
        final float constY = (a10*x) + (a11*y) + a12; //const part of Y coord
       	final float precalc2 = 2 * centerY;//const in inner loop quad. formula
        final float precalc3 =-2 * centerX;//const in inner loop quad. formula
        float X; // User space point X coordinate 
        float Y; // User space point Y coordinate
        float g;//value between 0 and 1 specifying position in the gradient
        float det; //determinant of quadratic formula (should always be >0)
        float currentToFocusSq;//sq distance from the current pt. to focus
        float intersectToFocusSq;//sq distance from the intersect pt. to focus
        float deltaXSq; //temp variable for a change in X squared.
        float deltaYSq; //temp variable for a change in Y squared.
        int indexer = off; //index variable for pixels array
        int i, j; //indexing variables for FOR loops
        int pixInc = w+adjust;//incremental index change for pixels array

        for (j = 0; j < h; j++) { //for every row
	    
            X = (a01*j) + constX; //constants from column to column
            Y = (a11*j) + constY;
	    
            //for every column (inner loop begins here)
            for (i = 0; i < w; i++) {	       			
	
                // special case to avoid divide by zero or very near zero
                if (((X-focusX)>-0.000001) &&
                    ((X-focusX)< 0.000001)) {		   
                    solutionX = focusX;
		    
                    solutionY = centerY;
		    
                    solutionY += (Y > focusY)?trivial:-trivial;
                }
		
                else {    
		    
                    //slope of the focus-current line
                    slope =   (Y - focusY) / (X - focusX);
		    
                    yintcpt = Y - (slope * X); //y-intercept of that same line
		    
                    //use the quadratic formula to calculate the intersection
                    //point		  
                    A = (slope * slope) + 1; 
		    
                    B =  precalc3 + (-2 * slope * (centerY - yintcpt));
		    
                    C =  constC + (yintcpt* (yintcpt - precalc2));
		    
                    det = (float)Math.sqrt((B * B) - ( 4 * A * C));
		    
                    solutionX = -B;
		    
                    //choose the positive or negative root depending
                    //on where the X coord lies with respect to the focus.
                    solutionX += (X < focusX)?-det:det;
		    
                    solutionX = solutionX / (2 * A);//divisor
		    
                    solutionY = (slope * solutionX) + yintcpt;
                }	                    	

                //calculate the square of the distance from the current point 
                //to the focus and the square of the distance from the 
                //intersection point to the focus. Want the squares so we can
                //do 1 square root after division instead of 2 before.

                deltaXSq = (float)solutionX - focusX;
                deltaXSq = deltaXSq * deltaXSq;

                deltaYSq = (float)solutionY - focusY;
                deltaYSq = deltaYSq * deltaYSq;

                intersectToFocusSq = deltaXSq + deltaYSq;

                deltaXSq = X - focusX;
                deltaXSq = deltaXSq * deltaXSq;

                deltaYSq = Y - focusY;
                deltaYSq = deltaYSq * deltaYSq;

                currentToFocusSq = deltaXSq + deltaYSq;

                //want the percentage (0-1) of the current point along the 
                //focus-circumference line
                g = (float)Math.sqrt(currentToFocusSq / intersectToFocusSq);

                //Get the color at this point
                pixels[indexer + i] = indexIntoGradientsArrays(g);
		
                X += a00; //incremental change in X, Y
                Y += a10;	
            } //end inner loop
            indexer += pixInc;
        } //end outer loop
    }


    /** Fill the raster, cycling the gradient colors when a point
     *  falls outside of the perimeter of the 100% stop circle. Use
     *  the anti-aliased gradient lookup.
     *
     *  This calculation first computes the intersection point of the line
     *  from the focus through the current point in the raster, and the
     *  perimeter of the gradient circle.
     * 
     *  Then it determines the percentage distance of the current point along
     *  that line (focus is 0%, perimeter is 100%). 
     *
     *  Equation of a circle centered at (a,b) with radius r:
     *  (x-a)^2 + (y-b)^2 = r^2
     *  Equation of a line with slope m and y-intercept b
     *  y = mx + b
     *  replacing y in the cirlce equation and solving using the quadratic
     *  formula produces the following set of equations.  Constant factors have
     *  been extracted out of the inner loop.
     * */
    private void antiAliasFillRaster(int pixels[], int off, 
                                     int adjust, 
                                     int x, int y, 
                                     int w, int h) {
        // Constant part of the C factor of the quadratic equation
        final double constC = 
            -(radiusSq) + (centerX * centerX) + (centerY * centerY);
        //coefficients of the quadratic equation (Ax^2 + Bx + C = 0)
       	final float precalc2 = 2 * centerY;//const in inner loop quad. formula
        final float precalc3 =-2 * centerX;//const in inner loop quad. formula

        //const part of X,Y coord (shifted to bottom left corner of pixel.
       	final float constX = (a00*(x-.5f)) + (a01*(y+.5f)) + a02;
        final float constY = (a10*(x-.5f)) + (a11*(y+.5f)) + a12;
        float X; // User space point X coordinate 
        float Y; // User space point Y coordinate
        int i, j; //indexing variables for FOR loops
        int indexer = off-1; //index variable for pixels array

        double [] prevGs = new double[w+1];
        double deltaXSq, deltaYSq;
        double solutionX, solutionY;
        double slope, yintcpt, A, B, C, det;
        double intersectToFocusSq, currentToFocusSq;
        double g00, g01, g10, g11;

        // Set X,Y to top left corner of first pixel of first row.
        X = constX - a01;
        Y = constY - a11;

        // Calc top row of g's.
        for (i=0; i <= w; i++) {
            // special case to avoid divide by zero or very near zero
            if (((X-focusX)>-0.000001) &&
                ((X-focusX)< 0.000001)) {		   
                solutionX = focusX;
                solutionY = centerY;
                solutionY += (Y > focusY)?trivial:-trivial;
            }
            else { 
                // Formula for Circle: (X-Xc)^2 + (Y-Yc)^2 - R^2 = 0
                // Formula line:        Y = Slope*x + Y0;
                // 
                // So you substitue line into Circle and apply
                // Quadradic formula.


                //slope of the focus-current line
                slope =   (Y - focusY) / (X - focusX);
		    
                yintcpt = Y - (slope * X); //y-intercept of that same line
		    
                //use the quadratic formula to calculate the intersection
                //point		  
                A = (slope * slope) + 1; 
		    
                B =  precalc3 + (-2 * slope * (centerY - yintcpt));
		    
                C =  constC + (yintcpt* (yintcpt - precalc2));
		    
                det = Math.sqrt((B * B) - ( 4 * A * C));
		    
                solutionX = -B;
		    
                //choose the positive or negative root depending
                //on where the X coord lies with respect to the focus.
                solutionX += (X < focusX)?-det:det;
		    
                solutionX = solutionX / (2 * A);//divisor
		    
                solutionY = (slope * solutionX) + yintcpt;
            }	                    	

            //calculate the square of the distance from the current point 
            //to the focus and the square of the distance from the 
            //intersection point to the focus. Want the squares so we can
            //do 1 square root after division instead of 2 before.
            deltaXSq = solutionX - focusX;
            deltaXSq = deltaXSq * deltaXSq;
        
            deltaYSq = solutionY - focusY;
            deltaYSq = deltaYSq * deltaYSq;
        
            intersectToFocusSq = deltaXSq + deltaYSq;
        
            deltaXSq = X - focusX;
            deltaXSq = deltaXSq * deltaXSq;
        
            deltaYSq = Y - focusY;
            deltaYSq = deltaYSq * deltaYSq;
        
            currentToFocusSq = deltaXSq + deltaYSq;
        
            //want the percentage (0-1) of the current point along the 
            //focus-circumference line
            prevGs[i] = Math.sqrt(currentToFocusSq / intersectToFocusSq);

            X += a00; //incremental change in X, Y
            Y += a10;	
        }

        for (j = 0; j < h; j++) { //for every row
	    
            // Set X,Y to bottom edge of pixel row.
            X = (a01*j) + constX; //constants from row to row
            Y = (a11*j) + constY;

            g10 = prevGs[0];
            // special case to avoid divide by zero or very near zero
            if (((X-focusX)>-0.000001) &&
                ((X-focusX)< 0.000001)) {		   
                solutionX = focusX;
                solutionY = centerY;
                solutionY += (Y > focusY)?trivial:-trivial;
            }
            else { 
                // Formula for Circle: (X-Xc)^2 + (Y-Yc)^2 - R^2 = 0
                // Formula line:        Y = Slope*x + Y0;
                // 
                // So you substitue line into Circle and apply
                // Quadradic formula.


                //slope of the focus-current line
                slope =   (Y - focusY) / (X - focusX);
		    
                yintcpt = Y - (slope * X); //y-intercept of that same line
		    
                //use the quadratic formula to calculate the intersection
                //point		  
                A = (slope * slope) + 1; 
		    
                B =  precalc3 + (-2 * slope * (centerY - yintcpt));
		    
                C =  constC + (yintcpt* (yintcpt - precalc2));
		    
                det = Math.sqrt((B * B) - ( 4 * A * C));
		    
                solutionX = -B;
		    
                //choose the positive or negative root depending
                //on where the X coord lies with respect to the focus.
                solutionX += (X < focusX)?-det:det;
		    
                solutionX = solutionX / (2 * A);//divisor
		    
                solutionY = (slope * solutionX) + yintcpt;
            }	                    	

            //calculate the square of the distance from the current point 
            //to the focus and the square of the distance from the 
            //intersection point to the focus. Want the squares so we can
            //do 1 square root after division instead of 2 before.
            deltaXSq = solutionX - focusX;
            deltaXSq = deltaXSq * deltaXSq;
        
            deltaYSq = solutionY - focusY;
            deltaYSq = deltaYSq * deltaYSq;
        
            intersectToFocusSq = deltaXSq + deltaYSq;
        
            deltaXSq = X - focusX;
            deltaXSq = deltaXSq * deltaXSq;
        
            deltaYSq = Y - focusY;
            deltaYSq = deltaYSq * deltaYSq;
                
            currentToFocusSq = deltaXSq + deltaYSq;
            g11 = Math.sqrt(currentToFocusSq / intersectToFocusSq);
            prevGs[0] = g11;
            
            X += a00; //incremental change in X, Y
            Y += a10;	

            //for every column (inner loop begins here)
            for (i=1; i <= w; i++) {
                g00 = g10;
                g01 = g11;
                g10 = prevGs[i];

                // special case to avoid divide by zero or very near zero
                if (((X-focusX)>-0.000001) &&
                    ((X-focusX)< 0.000001)) {		   
                    solutionX = focusX;
                    solutionY = centerY;
                    solutionY += (Y > focusY)?trivial:-trivial;
                }
                else { 
                    // Formula for Circle: (X-Xc)^2 + (Y-Yc)^2 - R^2 = 0
                    // Formula line:        Y = Slope*x + Y0;
                    // 
                    // So you substitue line into Circle and apply
                    // Quadradic formula.


                    //slope of the focus-current line
                    slope =   (Y - focusY) / (X - focusX);
		    
                    yintcpt = Y - (slope * X); //y-intercept of that same line
		    
                    //use the quadratic formula to calculate the intersection
                    //point		  
                    A = (slope * slope) + 1; 
		    
                    B =  precalc3 + (-2 * slope * (centerY - yintcpt));
		    
                    C =  constC + (yintcpt* (yintcpt - precalc2));
		    
                    det = Math.sqrt((B * B) - ( 4 * A * C));
		    
                    solutionX = -B;
		    
                    //choose the positive or negative root depending
                    //on where the X coord lies with respect to the focus.
                    solutionX += (X < focusX)?-det:det;
		    
                    solutionX = solutionX / (2 * A);//divisor
		    
                    solutionY = (slope * solutionX) + yintcpt;
                }	                    	

                //calculate the square of the distance from the current point 
                //to the focus and the square of the distance from the 
                //intersection point to the focus. Want the squares so we can
                //do 1 square root after division instead of 2 before.
                deltaXSq = solutionX - focusX;
                deltaXSq = deltaXSq * deltaXSq;
        
                deltaYSq = solutionY - focusY;
                deltaYSq = deltaYSq * deltaYSq;
        
                intersectToFocusSq = deltaXSq + deltaYSq;
        
                deltaXSq = X - focusX;
                deltaXSq = deltaXSq * deltaXSq;
        
                deltaYSq = Y - focusY;
                deltaYSq = deltaYSq * deltaYSq;
        
                currentToFocusSq = deltaXSq + deltaYSq;
                g11 = Math.sqrt(currentToFocusSq / intersectToFocusSq);
                prevGs[i] = g11;

                //Get the color at this point
                pixels[indexer+i] = indexGradientAntiAlias
                    ((float)((g00+g01+g10+g11)/4), 
                     (float)Math.max(Math.abs(g11-g00),
                                     Math.abs(g10-g01)));

                X += a00; //incremental change in X, Y
                Y += a10;	
            } //end inner loop
            indexer += (w+adjust);
        } //end outer loop
    }
}
