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_Cd, natom_Se;

//  Hardwired data for now
    double aa = 6.08; // CdSe

    double a_radius = 7.5, b_radius = 7.5, c_radius = 7.5;
    double exponent = 2.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, Cd_atoms, Se_atoms, Cd_in_basis, Se_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);

    Cd_atoms = new double [maxdim][3];
    Se_atoms = new double [maxdim][3];
    Cd_in_basis = new double [1][3];
    Se_in_basis = new double [1][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
   Cd_in_basis[0][0] = 0.0000 ; Cd_in_basis[0][1] = 0.000  ; Cd_in_basis[0][2] = 0.000;
   Se_in_basis[0][0] = 0.25*aa; Se_in_basis[0][1] = 0.25*aa; Se_in_basis[0][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_Cd = -1;
      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] = Cd_in_basis[0][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_Cd += 1;
              if (natom_Cd < maxdim) {
                Cd_atoms[natom_Cd][0] = tempo[0];
                Cd_atoms[natom_Cd][1] = tempo[1];
                Cd_atoms[natom_Cd][2] = tempo[2];

              }
            }
          }
        }
      }
  
      natom_Se = -1;
      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] = Se_in_basis[0][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_Se += 1;
              if (natom_Se < maxdim) {
                Se_atoms[natom_Se][0] = tempo[0];
                Se_atoms[natom_Se][1] = tempo[1];
                Se_atoms[natom_Se][2] = tempo[2];

              }
            }
          }
        }
      }
  
      p.printf ("%4d\n\n", natom_Cd + natom_Se);
      Data_out.writeInt(   natom_Cd + natom_Se);
  
      for( i=0; i < natom_Cd; i++) {
        p.printf ("Cd  %.4f  %.4f  %.4f\n",
                   Cd_atoms[i][0], Cd_atoms[i][1],Cd_atoms[i][2]);
      
        Data_out.writeUTF("Cd");
        Data_out.writeDouble(Cd_atoms[i][0]);
        Data_out.writeDouble(Cd_atoms[i][1]);
        Data_out.writeDouble(Cd_atoms[i][2]);
      }
      for( i=0; i < natom_Se; i++) {
        p.printf ("Se  %.4f  %.4f  %.4f\n",
                   Se_atoms[i][0], Se_atoms[i][1],Se_atoms[i][2]);
      
        Data_out.writeUTF("Se");
        Data_out.writeDouble(Se_atoms[i][0]);
        Data_out.writeDouble(Se_atoms[i][1]);
        Data_out.writeDouble(Se_atoms[i][2]);
      }
           
      p.close();
      Data_out.close();
    }
    catch (Exception e)
    {
      System.err.println ("Error writing to file");
    }


  }
}