//***********************************************************************************
// Generalized Julia set by complex iteration: f(z)=z^k+c (k Real; c Complex user supplied_
//
// Given complex value c and escape radius R (satisfying R(R-1)>=|c|) the Julia set J(c) is
// formed by the set of all Complex numbers for which
//
// J(c) = { z : for all n, z_{n+1}=(z_n)^k+c satisfies |z_{n+1}|<=R with initial value z_0=z}
//
// If z is not in J(c) then there will be some Whole number k with |z_k|>R (satarting z_0=z).
// This value can be used to determine a colour for the Complex number z in the Complex plane, all
// z in J(c) being coloured black (RGB=000).
//***************************************************************************
import java.awt.*;
import java.awt.event.*;
import java.awt.image.*;
import java.io.*;
import javax.imageio.*;
import javax.swing.*;
import java.util.*;
public class GenJulia
  {
  private static Scanner input = new Scanner(System.in);
  public static void main(String[] args) throws IOException
    {
    int ht = 10000;
    int wd = 10000;
    double R=0.0;  // Escape region
    System.out.println("Positive Real power k: for z^k+c");
    double ksize = input.nextDouble();
    System.out.println("Real c value: (for z_{n+1} = (z_n)^k+c)");
    double rec = input.nextDouble();
    System.out.println("Im c value: (for z_{n+1} = (z_n)^k+c)");
    double imc = input.nextDouble();
    double csize = Math.sqrt(rec*rec+imc*imc);
    System.out.println("Escape radius R (should satisfy R(R-1)>=|c|):");
    R = input.nextDouble();
    if (ksize<0) { R=2; } else {
    while (((Math.pow(R,ksize)-R)<csize))
      {
      System.out.println("Escape radius R (should satisfy R^k-R>=|c|):");
      R = input.nextDouble();
      };
    };
    double re_low = (-R);
    double re_high = R;
    double im_low = (-R);
    double im_high = R;
    if (R<2)
      {
       re_low = (-2.5);
       re_high = 3.5;
       im_low = (-2.5);
       im_high = 2.5;
      };
    double rez;
    double imz;
    double trez;
    double timz;
    double xnew;
    double ynew;
    int drop_out = 100;   // if |z_n|<=R for all n<=1000 assume z is in J(c)
    int tests_made;
    double re_inc = (re_high-re_low)/((double)wd);
    double im_inc = (im_high-im_low)/((double)ht);
    int xc;
    int yc;
    int colours = 128;
    int[] seen = new int[drop_out+1];
    int[][] colouring = new int[ht][wd];
    int[] palette = new int[colours];
    boolean inside;
    File CompImage = new File("GenJuliaSet"+"k"+ksize+"Re"+rec+"Im"+imc+".jpg");
    BufferedImage img2 = new BufferedImage(wd,ht, BufferedImage.TYPE_INT_RGB );
    // Set (default) colour palette
    // Complex values assigned one of 512 colours depending on exit time: resident (ie in J(c)) coloured Black;
    palette[0]=(255);
    for (int clr=1; clr<colours; clr++)
      {
      // To experiment with different colour templates, modify the 3 lines below:
      // (value)<<16 is the RED component in RGB
      // (value)<<8  is the GREEN
      // (value)     is the BLUE
      // Different functions of the index can be used, but is necessary to ensure component value is 8 bits, hence
      // modulo 256 operation
      int REDQ = ((clr+1)*(clr)*(clr-1))%256;  // How much Red
      int GREENQ = ((clr)*(clr))%256;  // How much Green
      int BLUEQ =((clr)*(clr))%256;  // How much Blue
         palette[clr] = (REDQ<<16)|(GREENQ<<8)|(BLUEQ);
      };
      for (int ctr=0; ctr<drop_out+1; ctr++)
        {
        seen[ctr]=0;
        };
    //************************************************************
    // Main process: iterate z_{n+1} = z_n^2+(rc,ic)
    // Classify (x,y) in Complex plane by how long to exit
    //************************************************************
    for (xc=0; xc<wd; xc++)
      {
      for (yc=0; yc<ht; yc++)
        {
        // Next z moving through Complex Plane
        rez = re_low+(double)(xc)*re_inc;
        imz = im_low+(double)(yc)*im_inc;
        tests_made=0;
        // Initiate z_0
        rez=rez+rec; imz=imz+imc;
        inside = (((rez*rez+imz*imz)<R*R));
        // Iterate z_{n+1}=(z_n)^2+c until either escape radius exceeded or z in J(c)
        while ((tests_made<drop_out)&&(inside))  
          {
          xnew = Math.pow(rez*rez+imz*imz,ksize/2.0)*Math.cos(ksize*Math.atan2(imz,rez))+rec;
          ynew = Math.pow(rez*rez+imz*imz,ksize/2.0)*Math.sin(ksize*Math.atan2(imz,rez))+imc;
          rez = xnew ; imz=ynew;
          inside = (((rez*rez+imz*imz)<R*R));
          tests_made++;
          };
          colouring[xc][yc]=tests_made;
          seen[tests_made]++;
        };
       };
        // Now map from palette colours to escape times 
        int next_col=0;
        for (int map=0; map<drop_out+1; map++)
          {
          if (seen[map]>0)
            {
            seen[map]=next_col; next_col=(next_col+1)%(colours);
            };
          };
          // Finally generate Image
          for (xc=0; xc<wd; xc++)
            {
            for (yc=0; yc<ht; yc++)
              {
              if (colouring[xc][yc]>=drop_out)
                {
                img2.setRGB(xc, ht-1-yc, 0); // Inside Julia set: colour BLACK
                }
             else
               {
               int coll = colouring[xc][yc];
               img2.setRGB(xc, ht-1-yc,palette[seen[coll]]); // Outside set: colour gradated on number needed
               };
            };
          };

   ImageIO.write(img2, "JPG", CompImage);
};
}