import java.io.*;
import java.lang.*;

public class ShapeDot2jmol {
  public static void main( String[] args ) {
    int i,j,k, icoord, iatom, iatom_in_cell;
    int natom;

//  Hardwired data for now
    double aa = 3.567; // Carbon

    double a_radius = 3.8, b_radius = 3.8, c_radius = 3.8;
    double exponent = 3.0;
    double anglex = 0.0; // Math.PI/4.0; 
    double angley = 0.0; // Math.PI/4.0;
    double anglez = 0.0; // Math.PI/4.0;

    double nn;

    double  Xcell, Ycell, Zcell; 
    double radius;

    double x,y,z;
    double [][] basis, all_atoms, atoms_in_basis; 
    double []   tempo, tempox, tempoy, tempoz;
    FileOutputStream out, FileData_out; // declare a file output object
    FileOutputStream parsec_out; // declare a file output object
    PrintStream p; // declare a print stream object

//  simple clear screen
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
    System.out.print("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");

    System.out.println("This program generates the core structure of quantum dot");
    System.out.println("with superellipsoid shape (Lamé curve)");
    System.out.println("|x/a|^n + |y/b|^n +|z/c|^n <= 1");
    System.out.println("Written by carrier@cs.umn.edu\n");

// java.lang.Math.pow(a,b);

    Xcell = 4.0*Math.ceil(a_radius/aa);
    Ycell = 4.0*Math.ceil(b_radius/aa);
    Zcell = 4.0*Math.ceil(c_radius/aa);

    int NXcell = (int) Xcell;
    int NYcell = (int) Ycell;
    int NZcell = (int) Zcell;

    int maxdim = (NXcell)*(NYcell)*(NZcell);
    System.out.printf("%4d X %4d X %4d = %4d atoms max\n",NXcell, NYcell, NZcell, maxdim);

    all_atoms = new double [maxdim][3];
    atoms_in_basis = new double [2][3];

    basis    = new double [3][3];
    tempo   = new double [3];
    tempox  = new double [3];
    tempoy  = new double [3];
    tempoz  = new double [3];

// The basis
   basis[0][0] = 0.0000; basis[0][1] = aa/2.0; basis[0][2] = aa/2.0;
   basis[1][0] = aa/2.0; basis[1][1] = 0.0000; basis[1][2] = aa/2.0;
   basis[2][0] = aa/2.0; basis[2][1] = aa/2.0; basis[2][2] = 0.0000;
    
// The list of atoms in the unit cell
   atoms_in_basis[0][0] = 0.0000; atoms_in_basis[0][1] = 0.000; atoms_in_basis[0][2] = 0.000;
   atoms_in_basis[1][0] = 0.25*aa; atoms_in_basis[1][1] = 0.25*aa; atoms_in_basis[1][2] = 0.25*aa;
    

// print the basis
    System.out.println("The lattice basis");
    System.out.printf( "%2.12f %2.12f %2.12f\n",basis[0][0], basis[0][1], basis[0][2] );
    System.out.printf( "%2.12f %2.12f %2.12f\n",basis[1][0], basis[1][1], basis[1][2] );
    System.out.printf( "%2.12f %2.12f %2.12f\n",basis[2][0], basis[2][1], basis[2][2] );

// prints to file out
    try
    {
      out = new FileOutputStream("quantumdot.xyz"); 
      p = new PrintStream( out ); 

      FileData_out = new FileOutputStream("quantumdot.dat");
      DataOutputStream Data_out = new DataOutputStream( FileData_out );

      natom = -1;
            for (iatom_in_cell=0; iatom_in_cell <= 1; iatom_in_cell++){
      for( i=-NXcell ; i<=NXcell; i++ ) {
        for( j=-NYcell ; j<=NYcell; j++ ) {
          for( k=-NZcell ; k<=NZcell; k++ ) {
              for(icoord=0 ; icoord<=2; icoord++ ) {
                tempo[icoord] = atoms_in_basis[iatom_in_cell][icoord] +
                         i*basis[0][icoord] + j*basis[1][icoord]+ k*basis[2][icoord];
              }
              // rotation Rx(alpha)
              tempox[0] =  tempo[0];
              tempox[1] =  Math.cos(anglex)*tempo[1] + Math.sin(anglex)*tempo[2];
              tempox[2] = -Math.sin(anglex)*tempo[1] + Math.cos(anglex)*tempo[2];

              // rotation Ry(alpha)
              tempoy[0] =  Math.cos(angley)*tempox[0] - Math.sin(angley)*tempox[2];
              tempoy[1] = tempox[1];
              tempoy[2] =  Math.sin(angley)*tempox[0] + Math.cos(angley)*tempox[2];

              // rotation Rz(alpha)
              tempoz[0] =  Math.cos(anglez)*tempoy[0] + Math.sin(anglez)*tempoy[1];
              tempoz[1] = -Math.sin(anglez)*tempoy[0] + Math.cos(anglez)*tempoy[1];
              tempoz[2] = tempoy[2];

              tempo[0] = tempoz[0];
              tempo[1] = tempoz[1];
              tempo[2] = tempoz[2];

              radius = Math.pow(Math.abs(tempo[0]/a_radius),exponent) +
                       Math.pow(Math.abs(tempo[1]/b_radius),exponent) +
                       Math.pow(Math.abs(tempo[2]/c_radius),exponent);
              double one_exponent = 1.0 / exponent;
              radius = Math.pow(radius, one_exponent);

              if (radius <= 1.0) {
                natom += 1;
                if (natom < maxdim) {
                  all_atoms[natom][0] = tempo[0];
                  all_atoms[natom][1] = tempo[1];
                  all_atoms[natom][2] = tempo[2];

                }
              }
            }
          }
        }
      }

      p.printf ("%4d\n\n", natom);
      Data_out.writeInt(natom);

      for( i=0; i < natom; i++) {
        p.printf ("C   %.4f  %.4f  %.4f\n",
                   all_atoms[i][0], all_atoms[i][1],all_atoms[i][2]);
        Data_out.writeUTF("C ");
        Data_out.writeDouble(all_atoms[i][0]);
        Data_out.writeDouble(all_atoms[i][1]);
        Data_out.writeDouble(all_atoms[i][2]);
      }
         
      p.close();
      Data_out.close();
    }
    catch (Exception e)
    {
      System.err.println ("Error writing to file");
    }


  }
}