#include <iostream>
using namespace std;

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <qfile.h>
#include <qtextstream.h>
#include <qstring.h>
#include <qvaluevector.h>

#define DIPOLAR_SUCCESS                  1
#define DIPOLAR_FILE_NOT_FOUND          -1
#define DIPOLAR_FILE_OPEN_FAILED        -2
#define DIPOLAR_WRONG_INPUT_FILE_FORMAT -3

#define CALC_POINT  0
#define CALC_PLANE  1
#define CALC_SPHERE 2

// The half the number of unit cells used to calculate the dipolar field.
// This should be absorbed into the input file in the future.
#define DIPOLAR_NO_UNIT_CELLS  50

#define PI 3.141592653589793238462643383279502884197

// where is this number from? Check units and comment!
#define UNITS 9.2741

//-----------------------------------------------------------------------
/**
 * <p>
 */
typedef struct {
  int     line_no;
  QString data;
} input_data;

//-----------------------------------------------------------------------
/**
 * <p>
 */
typedef struct {
  double pos[3];
  int    type;
} atom_pos;

//-----------------------------------------------------------------------
/**
 * <p>
 */
typedef struct {
  double component[3];
} mag_moment;

//-----------------------------------------------------------------------
/**
 * <p>
 */
typedef struct {
  int type;
  double dipoleMatrix[3][3];
  double field[3];
} dipole_type_info;

//-----------------------------------------------------------------------
/**
 * <p>
 * 
 * <p>calcTag: 0=single muon position, 1=dipole field in a plain, 2=dipole field on a sphere
 */
typedef struct {
  int calcTag; /// tag showing what to be calculated
  int noAtoms;
  int noMagTypes;
  QValueVector<atom_pos> atom;
  double latticeParam[3];
  double angle[3];
  double lorentzSphereRadius;
  QValueVector<mag_moment> magMoment;
  double muPos[3];
  double planeFix[3];    /// vector fixing a plane in the crystal coordinate system
  double planeNormal[3]; /// normal plane vector
  QValueVector<dipole_type_info> dipole;
} dipole_info;

//-----------------------------------------------------------------------
/**
 * <p>
 */
void dipolar_ordered_syntax()
{
  cout << endl << "syntax: dipolar_ordered input_file output_file [debug]";
  cout << endl << "  input_file: input file name";
  cout << endl << "    The structure of the input file is:";
  cout << endl << "     details to be written ...";
  cout << endl;
  cout << endl << " WARNING: all the coordinates (including the magnetic moments) needed";
  cout << endl << "          to be given in crystal structure coordinates. This is trivial";
  cout << endl << "          for structures with lattice angles all 90°, but might lead to";
  cout << endl << "          errors/confusions for lattice angles not equal to 90°"; 
  cout << endl << "     a typical example would be:";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "     %";
  cout << endl << "     % LSCO.inp";
  cout << endl << "     %";
  cout << endl << "     % input file to calculate magnetic dipol field distribution for an ordered state";
  cout << endl << "     %";
  cout << endl << "     % data from Japanese Journal of Applied Physics, Part 1 (1988) 27, 1132-1137  ";
  cout << endl << "     % ";
  cout << endl << "     % muon site from H.-U. Suter, Physica B 326 (2003) 329";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "     % magnetic lattice related input";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "                                    4  % number of magnetic atoms in the magnetic unit cell";
  cout << endl << "                                    2  % number of magnetic types in the magnetic unit cell";
  cout << endl << "";
  cout << endl << "         0.0000   0.0000   0.0000   1  % x, y, z, magnetic type";
  cout << endl << "         0.5000   0.5000   0.0000   2  % x, y, z, magnetic type";
  cout << endl << "         0.0000   0.5000   0.5000   1  % x, y, z, magnetic type";
  cout << endl << "         0.5000   0.0000   0.5000   2  % x, y, z, magnetic type";
  cout << endl << "";
  cout << endl << "         5.3535   5.4003  13.1441      % lattice parameters in angstroms (magnetic cell)";
  cout << endl << "        90.0000  90.0000  90.0000      % lattice angles in degrees";
  cout << endl << "";
  cout << endl << "      200.0000                         % Lorentz sphere radius in angstroms";
  cout << endl << "";
  cout << endl << "         0.6800   0.0000   0.0000      % magnetic moment for type 1";
  cout << endl << "        -0.6800   0.0000   0.0000      % magnetic moment for type 2";
  cout << endl << "";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "     % what to be calculated specifically follows now ";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "";
  cout << endl << "         point                         % tag showing that a single muon position to be calculated";
  cout << endl << "         0.1478   0.0000   0.1846      % muon site in lattice coordinates";
  cout << endl << "";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "     % end";
  cout << endl << "     %--------------------------------------------------------------------------------";
  cout << endl << "  output_file: output file name";
  cout << endl << endl;
}

//-----------------------------------------------------------------------
/**
 * <p>
 */
int dipolar_ordered_read_input_file(char *file_name, dipole_info *info)
{
  QFile      file(file_name);
  input_data line_data;
  QValueVector<input_data> data;
  
  bool       error = false;
  char       line[256];
  QString    qline;
  int        status;
  int        param_switch = 0;
  int        line_no = 0;
  atom_pos   atom;
  mag_moment magMoment;

  // check if the file exists
  if (!file.exists())
    return DIPOLAR_FILE_NOT_FOUND;

  // open file
  if (!file.open(IO_ReadOnly))
    return DIPOLAR_FILE_OPEN_FAILED;

  // read all the data
  int comment_pos;  
  while (file.readLine(&line[0], sizeof(line)) && !file.atEnd()) {
    qline = line;
    line_no++;
    // ignore empty and comment lines
    if (qline.stripWhiteSpace().isEmpty() || qline.startsWith("%")) 
      continue;
    // strip comments within data line
    comment_pos = qline.find("%");
    if (comment_pos>0)
      qline = qline.remove(comment_pos-1, qline.length()-comment_pos);
    qline = qline.lower();
    line_data.line_no = line_no;
    line_data.data = qline;
    data.append(line_data);
  }
  
  file.close();  

  // analyze data
  QValueVector<input_data>::Iterator it;
  for (it = data.begin(); it != data.end(); ++it) {
    switch (param_switch) {
      case 0: // no of magnetic atoms in the magnetic unit cell
        status = sscanf(it->data.latin1(), "%d", &info->noAtoms);
        if (status != 1) // likely to have wrong input file format
          error = true;
        param_switch++;
        break;
      case 1: // no of magnetic types in the magnetic unit cell
        status = sscanf(it->data.latin1(), "%d", &info->noMagTypes);
        if (status != 1) // likely to have wrong input file format
          error = true;
        param_switch++;
        break;
      case 2: // get all the atomic positions and the associated magnetic type
        status = sscanf(it->data.latin1(), "%lf %lf %lf %d", 
                        &atom.pos[0], &atom.pos[1], &atom.pos[2], &atom.type);
        if (status != 4) { // likely to have wrong input file format
          error = true;
        } else {
          info->atom.append(atom);
        }
        // check if all atom position are read
        if (info->atom.count() == (unsigned int)info->noAtoms) {
          param_switch++;
        }
        break;
      case 3: // lattice parameters
        status = sscanf(it->data.latin1(), "%lf %lf %lf", 
                        &info->latticeParam[0], &info->latticeParam[1], &info->latticeParam[2]);
        if (status != 3) // likely to have wrong input file format
          error = true;
        param_switch++;
        break;
      case 4: // lattice angles
        status = sscanf(it->data.latin1(), "%lf %lf %lf", 
                   &info->angle[0], &info->angle[1], &info->angle[2]);
        if (status != 3) // likely to have wrong input file format
          error = true;
        param_switch++;
        break;
      case 5: // Lorentz sphere radius in (A)
        status = sscanf(it->data.latin1(), "%lf", &info->lorentzSphereRadius);
        if (status != 1) // likely to have wrong input file format
          error = true;
        param_switch++;
        break;
      case 6: // get all components of the magnetic types
        status = sscanf(it->data.latin1(), "%lf %lf %lf", 
                    &magMoment.component[0], &magMoment.component[1], &magMoment.component[2]);
        if (status != 3) { // likely to have wrong input file format
          error = true;
        } else {
          info->magMoment.append(magMoment);
        }
        // check if all magnetic moments are read
        if (info->magMoment.count() == (unsigned int)info->noMagTypes) {
          param_switch++;
        }
        break;
      case 7: // tag what to be calculated
        if (it->data.contains("point") > 0) {
          info->calcTag = CALC_POINT;
        } else if (it->data.contains("plane") > 0) {
          info->calcTag = CALC_PLANE;
        } else if (it->data.contains("sphere") > 0) {
          info->calcTag = CALC_SPHERE;
        } else {
          error = true;
        }
        param_switch++;
        break;
      case 8: // depends on the calculation tag
        switch (info->calcTag) {
          case CALC_POINT: // get the muon position
            status = sscanf(it->data.latin1(), "%lf %lf %lf", 
                      &info->muPos[0], &info->muPos[1], &info->muPos[2]);
            if (status != 3)
              error = true;
            break;
          case CALC_PLANE:
            status = sscanf(it->data.latin1(), "%lf %lf %lf",
                      &info->planeFix[0], &info->planeFix[1], &info->planeFix[2]);
            if (status != 3)
              error = true;
            break;
          case CALC_SPHERE:
            break;
          default:
            break;            
        }
        param_switch++;
        break;
      case 9: // depends on the calculation tag
        switch (info->calcTag) {
          case CALC_PLANE:
            status = sscanf(it->data.latin1(), "%lf %lf %lf",
                      &info->planeNormal[0], &info->planeNormal[1], &info->planeNormal[2]);
            if (status != 3)
              error = true;          
            break;
          default:
            break;
        }
        param_switch++;
        break;      
      default:
        break;
    }
    
    if (error) {
      cout << endl << "ERROR in line " << it->line_no;
      return DIPOLAR_WRONG_INPUT_FILE_FORMAT;
    }
  }

  return DIPOLAR_SUCCESS;
}

//-----------------------------------------------------------------------
/**
 * <p>
 */
void dipolar_ordered_point_calc(dipole_info *info)
{
  double distance_vector[3];
  double direction_vector[3];
  double distance;  ///< distance from an atom to the muon
  double distance2; ///< squared distance form an atom to the muon
  double idist3;    ///< 1/distance^3
  double dcos[3];   ///< directional cosines
  
  // init dipole vector structure
  dipole_type_info dipole;
  for (int ii=0; ii<3; ii++) {
    for (int jj=0; jj<3; jj++) {
      dipole.dipoleMatrix[ii][jj] = 0.0; // init dipole tensor
    } 
    dipole.field[ii] = 0.0; // init dipole field
  }
  for (int i=0; i<info->noMagTypes; i++) {
    dipole.type = i+1; // init magnetic type
    info->dipole.append(dipole); // append to the vector
  } 
  
  // calculate the directional cosines
  for (int i=0; i<3; i++) {
    dcos[i] = cos(info->angle[i]*PI/180.0);
  }
  
  // triple loop over the x, y, z unit cells
  for (int uc_x = -DIPOLAR_NO_UNIT_CELLS; uc_x <= DIPOLAR_NO_UNIT_CELLS; uc_x++) {
    for (int uc_y = -DIPOLAR_NO_UNIT_CELLS; uc_y <= DIPOLAR_NO_UNIT_CELLS; uc_y++) {
      for (int uc_z = -DIPOLAR_NO_UNIT_CELLS; uc_z <= DIPOLAR_NO_UNIT_CELLS; uc_z++) {
        // iterate all atoms in the unit cell
        QValueVector<atom_pos>::iterator it;
        for (it = info->atom.begin(); it != info->atom.end(); ++it) {
        
          // calculate the distance vector of the magnetic atom to the muon
          distance_vector[0] = info->latticeParam[0]*((it->pos[0]+uc_x) - info->muPos[0]); // "x" component
          distance_vector[1] = info->latticeParam[1]*((it->pos[1]+uc_y) - info->muPos[1]); // "y" component
          distance_vector[2] = info->latticeParam[2]*((it->pos[2]+uc_z) - info->muPos[2]); // "z" component
          
          // calculate the real distance between the atom and the muon
          distance2 = distance_vector[0]*distance_vector[0]+
                      distance_vector[1]*distance_vector[1]+
                      distance_vector[2]*distance_vector[2]+
                      distance_vector[0]*distance_vector[1]*dcos[0]+
                      distance_vector[0]*distance_vector[2]*dcos[1]+
                      distance_vector[1]*distance_vector[2]*dcos[2];
          distance = sqrt(distance2);
          
          // calculate the direction vector
          direction_vector[0] = distance_vector[0]/distance;
          direction_vector[1] = distance_vector[1]/distance;
          direction_vector[2] = distance_vector[2]/distance;
                     
          // check if it is within the Lorentz sphere, if yes calculate the dipole field tensor
          if (info->lorentzSphereRadius - distance > 0.0) {
            idist3 = 1.0/(distance*distance2);
            for (int ii=0; ii<3; ii++) {
              for (int jj=0; jj<3; jj++) {
                if (ii==jj) { // diagonal elements
                  info->dipole.at(it->type-1).dipoleMatrix[ii][jj] += idist3*(3*direction_vector[ii]*direction_vector[ii]-1); 
                } else { // off diagonal elements
                  info->dipole.at(it->type-1).dipoleMatrix[ii][jj] += idist3*(3*direction_vector[ii]*direction_vector[jj]);
                }
              }         
            }         
          }
          
        }
      }
    }
  }

  // calculate the field
  for (int i=0; i<info->noMagTypes; i++) { // calculation for each magnetic type in the lattice
    for (int k=0; k<3; k++) {
      for (int l=0; l<3; l++) {
        info->dipole.at(i).field[k] += info->dipole.at(i).dipoleMatrix[k][l]*info->magMoment.at(i).component[l];
      }
    }
  }
}

//-----------------------------------------------------------------------
/**
 * <p>
 */
void dipolar_ordered_plane_calc(dipole_info *info)
{
  cout << endl << "SORRY: not yet implemented ..." << endl;
}

//-----------------------------------------------------------------------
/**
 * <p>
 */
void dipolar_ordered_sphere_calc(dipole_info *info)
{
  cout << endl << "SORRY: not yet implemented ..." << endl;
}

//-----------------------------------------------------------------------
/**
 * <p>
 */
int dipolar_ordered_write_output_file(char *input_file_name, char *output_file_name, dipole_info *info)
{
  char str[128];
  // copy the input file to the output file
  sprintf(str, "cp %s %s", input_file_name, output_file_name);
  system(str);

  // open the new output file and append the results
  QFile file(output_file_name);

  if (!file.open(IO_WriteOnly | IO_Append))
    return DIPOLAR_FILE_OPEN_FAILED;

  QTextStream fout(&file);

  double field_tot[3];
  for (int i=0; i<3; i++)
    field_tot[i] = 0.0;
    
  fout << endl << "% --------- begin output data -----------------";
  fout << endl;
  fout << endl << " ALL THE COMPONENTS ARE GIVEN AS COMPONENTS OF THE CRYSTAL, i.e.";
  fout << endl << " vec_B = B_a vec_a + B_b vec_b + B_c vec_c, etc., where e.g. vec_a is the";
  fout << endl << " normalized vector along the a-axis of the crystal etc.";
  fout << endl << " This needs to be kept in mind, especially for triklin, monoklin, ..., structures";
  fout << endl;
  QValueVector<dipole_type_info>::iterator it;
  int i=1;
  for (it = info->dipole.begin(); it != info->dipole.end(); ++it) {
    fout << endl;
    fout << endl << " THE DIPOLAR MATRIX FOR THE TYPE " << i++ << " in (kG/mu_B)";
    fout << endl;
    fout << endl << "  D |      1         2         3";
    fout << endl << "  --|-----------------------------------";
    for (int j=0; j<3; j++) {
      sprintf(str, "  %d |  %+2.4f   %+2.4f   %+2.4f  ", 
              j+1, UNITS*it->dipoleMatrix[j][0], UNITS*it->dipoleMatrix[j][1], UNITS*it->dipoleMatrix[j][2]);
      fout << endl << str;
      field_tot[j] += UNITS*it->field[j];
    }
    fout << endl;
    fout << endl << " trace = " << UNITS*(it->dipoleMatrix[0][0]+it->dipoleMatrix[1][1]+it->dipoleMatrix[2][2]);
    fout << endl;
    fout << endl << " the dipolar magnetic field generated by the type " << it->type << " ions is: ";
    fout << UNITS*sqrt(it->field[0]*it->field[0]+it->field[1]*it->field[1]+it->field[2]*it->field[2]) << " (kG)";
    fout << endl << " with the components: " << UNITS*it->field[0] << ", " << UNITS*it->field[1] << ", " << UNITS*it->field[2] << " (kG)";
    fout << endl;
    fout << endl << " -----------------------";
  }
  fout << endl;
  fout << endl << " THE TOTAL DIPOLAR MAGNETIC FIELD IS ";
  fout << sqrt(field_tot[0]*field_tot[0]+field_tot[1]*field_tot[1]+field_tot[2]*field_tot[2]) << " (kG)";
  fout << endl << " with the components: " << field_tot[0] << ", " << field_tot[1] << ", " << field_tot[2] << " (kG)";
  fout << endl;
  fout << endl << "% ---------  end output data  -----------------";

  file.close();

  return DIPOLAR_SUCCESS;
}

//-----------------------------------------------------------------------
/**
 * <p>
 */
int main(int argc, char *argv[])
{
  int debug = 0;
  dipole_info *info;
  int status;
  int i;

  // check syntax
  if (argc < 3) {
    dipolar_ordered_syntax();
    return 0;
  }

  if (argc == 4)
    sscanf(argv[3], "%d", &debug);

  // allocate memory for the info structure
  info = new dipole_info;

  // init info structure
  info->noAtoms = 0;
  info->noMagTypes = 0;
  info->latticeParam[0] = 0.0;
  info->latticeParam[1] = 0.0;
  info->latticeParam[2] = 0.0;
  info->angle[0] = 0.0;
  info->angle[1] = 0.0;
  info->angle[2] = 0.0;
  info->muPos[0] = 0.0;
  info->muPos[1] = 0.0;
  info->muPos[2] = 0.0;
  info->lorentzSphereRadius = 0.0;

  // read input file and fill info structure
  status = dipolar_ordered_read_input_file(argv[1], info);
  if (status != DIPOLAR_SUCCESS) {
    switch (status) {
      case DIPOLAR_FILE_NOT_FOUND:
        cout << endl << "ERROR: couldn't find input file " << argv[1] << endl;
        break;
      case DIPOLAR_FILE_OPEN_FAILED:
        cout << endl << "ERROR: couldn't open input file " << argv[1] << endl;
        break;
      case DIPOLAR_WRONG_INPUT_FILE_FORMAT:
        cout << endl << "ERROR: wrong input file format ..." << endl;
        break;
      default:
        cout << endl << "unknown error ..." << endl;
        break;
    }
    dipolar_ordered_syntax();
    return 0;
  }

  if (debug) {
    cout << endl << "-----------------------";
    cout << endl << "[debug] no of magnetic atoms in the unit cell = " << info->noAtoms;
    cout << endl << "[debug] no of magnetic types in the unit cell = " << info->noMagTypes;
    cout << endl << "[debug] list of the atom position (A) and its type";
    cout << endl << "-----------------------";
    i = 1;
    for (QValueVector<atom_pos>::iterator it = info->atom.begin(); it != info->atom.end(); ++it) {
      cout << endl << "[debug] " << i++ << ": ";
      cout << it->pos[0] << ", ";
      cout << it->pos[1] << ", ";
      cout << it->pos[2] << ", ";
      cout << it->type;
    }
    cout << endl << "-----------------------";
    cout << endl << "[debug] lattice parameters (A) : ";
    cout << info->latticeParam[0] << ", ";
    cout << info->latticeParam[1] << ", ";
    cout << info->latticeParam[2];
    cout << endl << "[debug] angles (°): ";
    cout << info->angle[0] << ", ";
    cout << info->angle[1] << ", ";
    cout << info->angle[2];
    cout << endl << "[debug] Lorentz sphere radius (A): ";
    cout << info->lorentzSphereRadius;
    cout << endl << "[debug] magnetic moments of the various types";
    cout << endl << "-----------------------";
    i = 1;
    for (QValueVector<mag_moment>::iterator it = info->magMoment.begin(); it != info->magMoment.end(); ++it) {
      cout << endl << "[debug] " << i++ << ": ";
      cout << it->component[0] << ", ";
      cout << it->component[1] << ", ";
      cout << it->component[2] << ", ";
    }
    cout << endl << "-----------------------";
    if (info->calcTag == CALC_POINT) {
      cout << endl << "[debug] calculate the dipolar field in a single point of the lattice.";
      cout << endl << "[debug] muon position in lattice coordinates: ";
      cout << info->muPos[0] << ", ";
      cout << info->muPos[1] << ", ";
      cout << info->muPos[2];
    } else if (info->calcTag == CALC_PLANE) {
      cout << endl << "[debug] calculate the dipolar field in a plane of the lattice.";
      cout << endl << "[debug] vector to the plane: ";
      cout << info->planeFix[0] << ", ";
      cout << info->planeFix[1] << ", ";
      cout << info->planeFix[2];
      cout << endl << "[debug] plane normal vector: ";
      cout << info->planeNormal[0] << ", ";
      cout << info->planeNormal[1] << ", ";
      cout << info->planeNormal[2];
    } else if (info->calcTag == CALC_SPHERE) {
      cout << endl << "[debug] calculate the dipolar field on a sphere.";
    }
    cout << endl << "-----------------------";
    cout << endl << endl;
  }

  // do the calculus
  switch (info->calcTag) {
    case CALC_POINT:
      dipolar_ordered_point_calc(info);
      break;
    case CALC_PLANE:
      dipolar_ordered_plane_calc(info);
      break;
    case CALC_SPHERE:
      dipolar_ordered_sphere_calc(info);
      break;
    default:
      cout << endl << "SEVERE ERROR: you never should get to here!?!?!?! Sorry ;-(" << endl; 
      break;
  }

  // write the output file
  status = dipolar_ordered_write_output_file(argv[1], argv[2], info);
  if (status != DIPOLAR_SUCCESS) {
    cout << endl << "ERROR while trying to write to the output file " << argv[2];
    cout << endl << endl;
  }

  // clean up
  if (info) {
    info->atom.clear();
    delete info;
    info = 0;
  }

  return 1;
}