#include "dump2sq.h"

void SetMatrix(HEAD *head){
  double a, b, c;
  double alpha, beta, gamma;

  head->lmax = head->lx;
  if (head->lmax < head->ly) head->lmax = head->ly;
  if (head->lmax < head->lz) head->lmax = head->lz;
  head->lmin = head->lx;
  if (head->lmin > head->ly) head->lmin = head->ly;
  if (head->lmin > head->lz) head->lmin = head->lz;
  
  a = head->lx;
  b = head->ly;
  c = head->lz;

  alpha = head->alpha * M_PI/180.0;
  beta = head->beta * M_PI/180.0;
  gamma = head->gamma * M_PI/180.0;
  head->M[0][0] = a;
  head->M[0][1] = 0;
  head->M[0][2] = 0;
  head->M[1][0] = b*cos(gamma);
  head->M[1][1] = b*sin(gamma);
  head->M[1][2] = 0;
  head->M[2][0] = c*cos(beta);
  head->M[2][1] = c*(cos(alpha)-cos(beta)*cos(gamma))/sin(gamma);
  head->M[2][2] = c*sqrt(1+2*cos(alpha)*cos(beta)*cos(gamma)
		   - cos(alpha)*cos(alpha)
		   - cos(beta)*cos(beta)
		   - cos(gamma)*cos(gamma)) /sin(gamma);
  head->M[2][2] = c*sqrt(1+2*cos(alpha)*cos(beta)*cos(gamma)
                          - cos(alpha)*cos(alpha)
                          - cos(beta)*cos(beta)
                          - cos(gamma)*cos(gamma)) /sin(gamma);
  
  head->volume = head->M[0][0] * head->M[1][1] * head->M[2][2];
  head->rho = head->atoms / head->volume;
  
  /* 逆行列の計算 */
  double det_A;
  det_A = head->volume;
  head->M_[0][0] = (head->M[1][1]*head->M[2][2] - head->M[1][2]*head->M[2][1])/det_A;
  head->M_[1][0] = (head->M[1][2]*head->M[2][0] - head->M[1][0]*head->M[2][2])/det_A;
  head->M_[2][0] = (head->M[1][0]*head->M[2][1] - head->M[1][1]*head->M[2][0])/det_A;
  head->M_[0][1] = 0.0;
  head->M_[1][1] = (head->M[0][0]*head->M[2][2] - head->M[0][2]*head->M[2][0])/det_A;
  head->M_[2][1] = (head->M[0][1]*head->M[2][0] - head->M[0][0]*head->M[2][1])/det_A;
  head->M_[0][2] = 0.0;
  head->M_[1][2] = 0.0;
  head->M_[2][2] = (head->M[0][0]*head->M[1][1] - head->M[0][1]*head->M[1][0])/det_A;
}

void ToFractionalCoord(HEAD *head, DATA *data){
  int i;
  double x, y, z;

  for (i=0; i<head->atoms; i++){
    x = data[i].x;
    y = data[i].y;
    z = data[i].z;
    data[i].x = x*head->M_[0][0] + y*head->M_[1][0] + z*head->M_[2][0];
    data[i].y = x*head->M_[0][1] + y*head->M_[1][1] + z*head->M_[2][1];
    data[i].z = x*head->M_[0][2] + y*head->M_[1][2] + z*head->M_[2][2];
    data[i].x = data[i].x - floor(data[i].x);
    data[i].y = data[i].y - floor(data[i].y);
    data[i].z = data[i].z - floor(data[i].z);
  }
}

/* femteckのfort.10用 */
int ReadHeadXYZ(FILE *fp, HEAD *head){
  char buf[1024];
  int chk;

  /* ITEM: TIMESTEP */
  if(fgets(buf, 1024, fp)==NULL) return 1;
  chk = sscanf(buf, "%d", &head->atoms);
  if(fgets(buf, 1024, fp)==NULL) return 1;
  chk = sscanf(buf, "Step %d %lf %lf %lf %lf %lf %lf",
	       &head->step,
	       &head->lx, &head->ly, &head->lz,
	       &head->alpha, &head->beta, &head->gamma);
  if (chk == 4){
    head->alpha = 90.0;
    head->beta = 90.0;
    head->gamma = 90.0;
  }
  SetMatrix(head);
  return 0;
}

/* femteckのfort.10用 */
int ReadDataXYZ(FILE *fp, HEAD *head, DATA *data)
{
  int i;
  char buf[1024];
  
  for(i=0; i<head->atoms; i++){
    if(fgets(buf, sizeof(buf), fp)==NULL) return 1;
    sscanf(buf, "%s %lf %lf %lf",
	   data[i].elem, &data[i].x, &data[i].y, &data[i].z);
  }
  ToFractionalCoord(head, data);
  return 0;
}

int ReadHeadLammpstrj(FILE *fp, HEAD *head){
  char buf[1024];
  int chk;
  double a, b, c;
  double ly, lz;
  double xlo, xhi, xy;
  double ylo, yhi, xz;
  double zlo, zhi, yz;
  double x_low=0.0, x_high=0.0;
  double y_low=0.0, y_high=0.0;

  /* ITEM: TIMESTEP */
  if(fgets(buf, 1024, fp)==NULL) return 1;
  if(fgets(buf, 1024, fp)==NULL) return 1;
  /* ITEM: NUMBER OF ATOMS */
  if(fgets(buf, 1024, fp)==NULL) return 1;
  if(fgets(buf, 1024, fp)==NULL) return 1;
  chk = sscanf(buf, "%d", &head->atoms);
  if (chk != 1) return 1;
  /* ITEM: BOX BOUNDS pp pp pp */
  if(fgets(buf, 1024, fp)==NULL) return 1;
  if(fgets(buf, 1024, fp)==NULL) return 1;
  chk = sscanf(buf, "%lf %lf %lf", &xlo, &xhi, &xy);
  if(fgets(buf, 1024, fp)==NULL) return 1;  
  chk = sscanf(buf, "%lf %lf %lf", &ylo, &yhi, &xz);
  if(fgets(buf, 1024, fp)==NULL) return 1;    
  chk = sscanf(buf, "%lf %lf %lf", &zlo, &zhi, &yz);
  if (chk == 3){ /* triclinic */
    /* x */
    if (x_low> xy) x_low = xy;
    if (x_low > xz) x_low = xz;
    if (x_low > (xy+xz)) x_low = xy+xz;
    if (x_high < xy) x_high = xy;
    if (x_high < xz) x_high = xz;
    if (x_high < (xy+xz)) x_high = xy+xz;
    xlo = xlo - x_low;
    xhi = xhi - x_high;
    a = xhi - xlo;
    /* y */
    if (y_low > yz) y_low = yz;
    if (y_high < yz) y_high = yz;
    ylo = ylo - y_low;
    yhi = yhi - y_high;
    ly = yhi - ylo;
    b = sqrt(ly*ly + xy*xy);
    /* z */
    lz = zhi - zlo;
    c = sqrt(lz*lz + xz*xz + yz*yz);
    /* angle */
    head->alpha = acos((yz*ly+xy*xz)/ (b*c)) * 180.0/M_PI;
    head->beta = acos(xz/c) * 180.0/M_PI;
    head->gamma = acos(xy/b) * 180.0/M_PI;
    head->lx = a;
    head->ly = b;
    head->lz = c;
  }
  else{
    head->alpha = 90;
    head->beta = 90;
    head->gamma = 90;
    head->lx = xhi - xlo;
    head->ly = yhi - ylo;
    head->lz = zhi - zlo;
  }
  SetMatrix(head);
  return 0;
}

int ReadDataLammpstrj(FILE *fp, HEAD *head, DATA *data)
{
  int i;
  char buf[1024];

  if(fgets(buf, sizeof(buf), fp)==NULL) return 1;
  for(i=0; i<head->atoms; i++){
    if(fgets(buf, sizeof(buf), fp)==NULL) return 1;  
    sscanf(buf, "%d %d %s %lf %lf %lf",
           &data[i].id, &data[i].type,
           data[i].elem,
           &data[i].x, &data[i].y, &data[i].z);
  }
  ToFractionalCoord(head, data);
  return 0;
}