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

//  Hardwired data for now
    double aa = 6.124; // PbSe

// 7.7 gives a 3 bonds min structure
    double a_radius = 8.7, b_radius = 8.7, c_radius = 8.7;
    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, Pb_atoms, Se_atoms, Pb_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");

    Xcell = 5.0*Math.ceil(a_radius/aa);
    Ycell = 5.0*Math.ceil(b_radius/aa);
    Zcell = 5.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);

    Pb_atoms = new double [maxdim][3];
    Se_atoms = new double [maxdim][3];
    Pb_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
   Pb_in_basis[0][0] = 0.0000 ; Pb_in_basis[0][1] = 0.000  ; Pb_in_basis[0][2] = 0.000;
   Se_in_basis[0][0] = 0.50*aa; Se_in_basis[0][1] = 0.50*aa; Se_in_basis[0][2] = 0.50*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_Pb = -1;
      for( i=-NXcell-1 ; i<=NXcell+1; i++ ) {
        for( j=-NYcell-1 ; j<=NYcell+1; j++ ) {
          for( k=-NZcell-1 ; k<=NZcell+1; k++ ) {
            for(icoord=0 ; icoord<=2; icoord++ ) {
              tempo[icoord] = Pb_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_Pb += 1;
              if (natom_Pb < maxdim) {
//      p.printf ("Pb  %.4f  %.4f  %.4f\n", tempo[0], tempo[1],tempo[2]);
                Pb_atoms[natom_Pb][0] = tempo[0];
                Pb_atoms[natom_Pb][1] = tempo[1];
                Pb_atoms[natom_Pb][2] = tempo[2];
              }
            }
          }
        }
      }
  
      natom_Se = -1;
      for( i=-NXcell-1 ; i<=NXcell+1; i++ ) {
        for( j=-NYcell-1 ; j<=NYcell+1; j++ ) {
          for( k=-NZcell-1 ; k<=NZcell+1; 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) {
//      p.printf ("Se  %.4f  %.4f  %.4f\n", tempo[0], tempo[1],tempo[2]);
                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_Pb+1 + natom_Se+1);
      Data_out.writeInt(   natom_Pb+1 + natom_Se+1);
  
      for( i=0; i <= natom_Pb; i++) {
        p.printf ("Pb  %.4f  %.4f  %.4f\n",
                   Pb_atoms[i][0], Pb_atoms[i][1],Pb_atoms[i][2]);
      
        Data_out.writeUTF("Pb");
        Data_out.writeDouble(Pb_atoms[i][0]);
        Data_out.writeDouble(Pb_atoms[i][1]);
        Data_out.writeDouble(Pb_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");
    }


  }
}