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

# define PI 3.14159265358979323846264338328

typedef struct {
  double x1;
  double x2;
  double x3;
} POINT;

typedef struct {
  double a;
  double b;
  double c;
  double alpha;
  double beta;
  double gamma;
} CELL;         /* cell parameters */

typedef struct {
  double aa1;
  double aa2;
  double aa3;
  double bb1;
  double bb2;
  double bb3;
  double cc1;
  double cc2;
  double cc3;
} CELL_C;      /* coordinates of cell corners */

typedef struct {
  double phi;
  double chi;
  double psi;
  double cos_phi;
  double cos_chi;
  double cos_psi;
  double sin_phi;
  double sin_chi;
  double sin_psi;
} ROTATE;      /* rotation parameters */

typedef struct {
  int h;
  int k;
  int l;
} MILLER_INDEX;

POINT lattice_coord(CELL *, int);
CELL_C reciprocal_coord(CELL *);
void coordinate(CELL_C *, ROTATE *, MILLER_INDEX *, POINT *);
void rotateX(double, double, POINT *);
/* void rotateY(double, double, POINT *); */
void rotateZ(double, double, POINT *);
double distance2(POINT *);
POINT cross_prdct(POINT *, POINT *);
/* double dot_prdct(POINT *, POINT *); */
int do_cross(POINT *, double *, double *, double *, double *);
int check_range(double, double, int, int, double *, double);

int main(int argc, char **argv)
{
  CELL *cell;                         /* Cell parameter */
  CELL_C rcell;                       /* (Reciprocal) cell corner coordinate */
  double mosaic, c_mosaic, s_mosaic;  /* Mosaicity of sample crystal : isotropic : full width */
  double lambda;                      /* Wave length */
  double LMT_SPHR, resolim;           /* Resolution limit : flag_reso=0 -> resolim(degree)
                                                            flag_reso=1 -> resolim(x/lambda) , x <= 2.0 */
  ROTATE *PhiChiPsi;                  /* Rotation angle : Z-X-Z extrinsic rotation phi(Z)->chi(X)->psi(Z)
                                         Euler(alpha, beta, gamma) : alpha=psi, beta=chi, gamma=phi */
  MILLER_INDEX *MI;                   /* Miller index */
  int hlim, klim, llim;               /* Limits of Miller indices */
  POINT *Q;                    /* Coordinates */
  double R2, Scale, PX, PY, LMT2;
  FILE *diffractionRLP;
  int flag_reso, flag_mosaic;         /* flag */

  /* Initialize */
  cell = (CELL *)malloc(sizeof(CELL));
  Q = (POINT *)malloc(sizeof(POINT));
  PhiChiPsi = (ROTATE *)malloc(sizeof(ROTATE));
  MI = (MILLER_INDEX *)malloc(sizeof(MILLER_INDEX));

  --argc;
  ++argv;
  if ((argc == 7) || (argc == 10) || (argc == 12)|| (argc == 14)) {
    sscanf(*argv, "%lf", &cell->a); argv++;              /* a */
    sscanf(*argv, "%lf", &cell->b); argv++;              /* b */
    sscanf(*argv, "%lf", &cell->c); argv++;              /* c */
    cell->alpha = 90.0;                                      /* alpha */
    cell->beta = 90.0;                                      /* beta */
    cell->gamma = 90.0;                                      /* gamma */
    if ((argc == 10)||(argc == 12) || (argc == 14)) {
    sscanf(*argv, "%lf", &cell->alpha); argv++;              /* alpha */
    sscanf(*argv, "%lf", &cell->beta); argv++;              /* beta */
    sscanf(*argv, "%lf", &cell->gamma); argv++;              /* gamma */
      }
    sscanf(*argv, "%lf", &PhiChiPsi->phi); argv++;          /* phi */
    sscanf(*argv, "%lf", &PhiChiPsi->chi); argv++;          /* chi */
    sscanf(*argv, "%lf", &PhiChiPsi->psi); argv++;          /* psi */
    sscanf(*argv, "%lf", &lambda); argv++;                 /* lambda */
    flag_reso=0;
    resolim=45.0; /* Default : 2theta = 45.0 degree */
    if ((argc == 12) || (argc == 14)) {
      sscanf(*argv, "%d", &flag_reso); argv++;             /* flag_reso */
      sscanf(*argv, "%lf", &resolim); argv++;              /* resolution 2theta or |s|/lambda */ 
    }
    flag_mosaic=0;
    mosaic=1.0; /* Default : Mosaicity = 1.0 degree , Full width */
    if (argc == 14) {
      sscanf(*argv, "%d", &flag_mosaic); argv++;           /* flag_mosaic */
      sscanf(*argv, "%lf", &mosaic);                       /* mosaicity */
    }
  } else {
    fputs("Number of arguments is inadequate.\n", stderr);
    exit(1);
  }

  /* To check input parameters in the range */
  if (1 != check_range(0.0, 500.0, 0, 1, &cell->a, 1.0)) fputs("a is out of the range. Set to 1.0 .\n", stderr);
  if (1 != check_range(0.0, 500.0, 0, 1, &cell->b, 1.0)) fputs("b is out of the range. Set to 1.0 .\n", stderr);
  if (1 != check_range(0.0, 500.0, 0, 1, &cell->c, 1.0)) fputs("c is out of the range. Set to 1.0 .\n", stderr);
  if (1 != check_range(0.0, 180.0, 0, 0, &cell->alpha, 90.0)) fputs("Alpha is out of the range. Set to 90.0 .\n", stderr);
  if (1 != check_range(0.0, 180.0, 0, 0, &cell->beta, 90.0)) fputs("Beta is out of the range. Set to 90.0 .\n", stderr);
  if (1 != check_range(0.0, 180.0, 0, 0, &cell->gamma, 90.0)) fputs("Gamma is out of the range. Set to 90.0 .\n", stderr);
  if (1 != check_range(0.1, 10.0, 1, 1, &lambda, 1.0)) fputs("Lambda is out of the range. Set to 1.0 .\n", stderr);

  if (flag_reso == 1) {
    if (1 != check_range(0.0, 2.5, 0, 1, &resolim, 2.0)) fputs("Resolution is out of range. Set to 2.0 (1/lambda).\n", stderr);
    LMT_SPHR = resolim;
  } else {
    if (flag_reso == 0) {
      if (1 != check_range(0.0, 180.0, 0, 1, &resolim, 180.0)) fputs("Resolution is out of range. Set to 180.0 (degree).\n", stderr);
      LMT_SPHR = 2.0 * sin((resolim/2.0) * PI/180.0);
    } else {
      fputs("Resolution flag_reso is incorrect.\n", stderr);
      exit(1);
    }
  }

  if (flag_mosaic != 1) {
    if (flag_mosaic != 0) {
      fputs("Mosaicity flag is incorrect. Set to 1.0 degree.\n", stderr);
    }
    mosaic = 1.0;
    flag_mosaic = 0;
  } else {
    if (1 != check_range(0.0, 10.0, 0, 1, &mosaic, 1.0)) fputs("Mosaicity is out of the range. Set to 1.0\n", stderr);
  }


  Scale = 1.0 / lambda;
  LMT2 = LMT_SPHR * LMT_SPHR * Scale * Scale;
  /* Limits of Miller Indices : +20% margin */
  hlim = ceil(LMT_SPHR * cell->a * Scale * 1.2);
  klim = ceil(LMT_SPHR * cell->b * Scale * 1.2);
  llim = ceil(LMT_SPHR * cell->c * Scale * 1.2);
  /* Rotation Parameters */
  PhiChiPsi->cos_phi = cos(PhiChiPsi->phi * PI/180.0);
  PhiChiPsi->cos_chi = cos(PhiChiPsi->chi * PI/180.0);
  PhiChiPsi->cos_psi = cos(PhiChiPsi->psi * PI/180.0);
  PhiChiPsi->sin_phi = sin(PhiChiPsi->phi * PI/180.0);
  PhiChiPsi->sin_chi = sin(PhiChiPsi->chi * PI/180.0);
  PhiChiPsi->sin_psi = sin(PhiChiPsi->psi * PI/180.0);
  /* Conversion to Radian */
  cell->alpha *= PI/180.0;
  cell->beta *= PI/180.0;
  cell->gamma *= PI/180.0;
  mosaic *= PI/180.0;

  c_mosaic = cos(mosaic/2.0);
  s_mosaic = sin(mosaic/2.0);

  /* Calculation of reciprocal cell */
  rcell = reciprocal_coord(cell);

  diffractionRLP = fopen("diffractionRLP.dat","w");

  for (MI->h = -hlim; MI->h <= hlim; (MI->h)++) {
    for (MI->k = -klim; MI->k <= klim; (MI->k)++) {
      for (MI->l = -llim; MI->l <= llim; (MI->l)++) {
	coordinate(&rcell, PhiChiPsi, MI, Q);
	R2 = distance2(Q);
	if (R2 < LMT2) {
	  if (Q->x1 <= 0.0) {
	    if (do_cross(Q, &Scale, &R2, &c_mosaic, &s_mosaic)) {
	      if (Q->x1 > (- Scale)) {
		PX = - (Q->x2 / (Scale + Q->x1));
		PY = (Q->x3 / (Scale + Q->x1));
	      /* Diffraction occurs. On the plate camera. */
		fprintf(diffractionRLP, "%4d %4d %4d  %9.5f %9.5f %9.5f  %9.5f %9.5f %9.5f  1  1  %9.5f %9.5f %9.5f\n", MI->h, MI->k, MI->l, Q->x1, Q->x2, Q->x3, Scale, PX, PY, PhiChiPsi->phi, PhiChiPsi->chi, PhiChiPsi->psi);
	      } else {
	      /* Diffraction occurs. Not on the plate camera. */
		fprintf(diffractionRLP, "%4d %4d %4d  %9.5f %9.5f %9.5f  %9.5f 0.0 0.0  1  0  %9.5f %9.5f %9.5f\n", MI->h, MI->k, MI->l, Q->x1, Q->x2, Q->x3, Scale, PhiChiPsi->phi, PhiChiPsi->chi, PhiChiPsi->psi);
	      }
	    } else {
	      /* Diffraction does not occurs. */
		fprintf(diffractionRLP, "%4d %4d %4d  %9.5f %9.5f %9.5f  %9.5f 0.0 0.0  0  0  %9.5f %9.5f %9.5f\n", MI->h, MI->k, MI->l, Q->x1, Q->x2, Q->x3, Scale, PhiChiPsi->phi, PhiChiPsi->chi, PhiChiPsi->psi);
	    }
	  } else {
	      /* Diffraction does not occurs. */
		fprintf(diffractionRLP, "%4d %4d %4d  %9.5f %9.5f %9.5f  %9.5f 0.0 0.0  0  0  %9.5f %9.5f %9.5f\n", MI->h, MI->k, MI->l, Q->x1, Q->x2, Q->x3, Scale, PhiChiPsi->phi, PhiChiPsi->chi, PhiChiPsi->psi);
	  }
	}
      }
    }
  }
  fclose(diffractionRLP);
  return(0);
}

CELL_C reciprocal_coord(CELL *cell_ptr){
  CELL_C RCC;
  POINT Qa, Qb, Qc, RQa, RQb, RQc;
  double V;
  Qa = lattice_coord(cell_ptr, 1);
  Qb = lattice_coord(cell_ptr, 2);
  Qc = lattice_coord(cell_ptr, 3);
  V = cell_ptr->a * cell_ptr->b * cell_ptr->c * sqrt(1 - cos(cell_ptr->alpha)*cos(cell_ptr->alpha) - cos(cell_ptr->beta)*cos(cell_ptr->beta) - cos(cell_ptr->gamma)*cos(cell_ptr->gamma) * 2*cos(cell_ptr->alpha)*cos(cell_ptr->beta)*cos(cell_ptr->gamma));
  RQa = cross_prdct(&Qb, &Qc);
  RQb = cross_prdct(&Qc, &Qa);
  RQc = cross_prdct(&Qa, &Qb);
  RCC.aa1 = RQa.x1/V;
  RCC.aa2 = RQa.x2/V;
  RCC.aa3 = RQa.x3/V;
  RCC.bb1 = RQb.x1/V;
  RCC.bb2 = RQb.x2/V;
  RCC.bb3 = RQb.x3/V;
  RCC.cc1 = RQc.x1/V;
  RCC.cc2 = RQc.x2/V;
  RCC.cc3 = RQc.x3/V;
  return(RCC);
}

POINT lattice_coord(CELL *cell_ptr, int axis)
{
  POINT cart;
  double c, c1, c2, c3;
  if (axis == 1) {
    cart.x1 = cell_ptr->a;
    cart.x2 = 0.0;
    cart.x3 = 0.0;
  }
  if (axis == 2) {
    cart.x1 = cell_ptr->b * cos(cell_ptr->gamma);
    cart.x2 = cell_ptr->b * sin(cell_ptr->gamma);
    cart.x3 = 0.0;
  }
  if (axis == 3) {
    c = cell_ptr->c;
    c1 = c * cos(cell_ptr->beta);
    c2 = c * (cos(cell_ptr->alpha) - cos(cell_ptr->beta) * cos(cell_ptr->gamma)) / sin(cell_ptr->gamma);
    c3 = sqrt(c * c - c1 * c1 - c2 * c2);
    cart.x1 = c1;
    cart.x2 = c2;
    cart.x3 = c3;
  }
  return(cart);
}

void coordinate(CELL_C *RLc_ptr, ROTATE *rot, MILLER_INDEX *Index_ptr, POINT *Q_ptr)
{
  Q_ptr->x1 = RLc_ptr->aa1 * Index_ptr->h + RLc_ptr->bb1 * Index_ptr->k + RLc_ptr->cc1 * Index_ptr->l;
  Q_ptr->x2 = RLc_ptr->aa2 * Index_ptr->h + RLc_ptr->bb2 * Index_ptr->k + RLc_ptr->cc2 * Index_ptr->l;
  Q_ptr->x3 = RLc_ptr->aa3 * Index_ptr->h + RLc_ptr->bb3 * Index_ptr->k + RLc_ptr->cc3 * Index_ptr->l;
  rotateZ(rot->cos_phi, rot->sin_phi, Q_ptr);
  rotateX(rot->cos_chi, rot->sin_chi, Q_ptr);
  rotateZ(rot->cos_psi, rot->sin_psi, Q_ptr);
}

void rotateX(double cos_phi, double sin_phi, POINT *Q_ptr)
{
  double rx, ry, rz;
  rx = Q_ptr->x1;
  ry = cos_phi * (Q_ptr->x2) - sin_phi * (Q_ptr->x3);
  rz = sin_phi * (Q_ptr->x2) + cos_phi * (Q_ptr->x3);
  Q_ptr->x1 = rx;
  Q_ptr->x2 = ry;
  Q_ptr->x3 = rz;
}

/* void rotateY(double cos_chi, double sin_chi, POINT *Q_ptr)
{
  double rx, ry, rz;
  rx = cos_chi * (Q_ptr->x1) + sin_chi * (Q_ptr->x3);
  ry = Q_ptr->x2;
  rz = -sin_chi * (Q_ptr->x1) + cos_chi * (Q_ptr->x3);
  Q_ptr->x1 = rx;
  Q_ptr->x2 = ry;
  Q_ptr->x3 = rz;
  } */

void rotateZ(double cos_psi, double sin_psi, POINT *Q_ptr)
{
  double rx, ry, rz;
  rx = cos_psi * (Q_ptr->x1) - sin_psi * (Q_ptr->x2);
  ry = sin_psi * (Q_ptr->x1) + cos_psi * (Q_ptr->x2);
  rz = Q_ptr->x3;
  Q_ptr->x1 = rx;
  Q_ptr->x2 = ry;
  Q_ptr->x3 = rz;
}

double distance2(POINT *Q_ptr)
{
  double r2;
  r2 = Q_ptr->x1 * Q_ptr->x1 + Q_ptr->x2 * Q_ptr->x2 + Q_ptr->x3 * Q_ptr->x3;
  return(r2);
}

POINT cross_prdct(POINT *v1, POINT *v2)
{
POINT v3;
v3.x1 = v1->x2 * v2->x3 - v1->x3 * v2->x2;
v3.x2 = v1->x3 * v2->x1 - v1->x1 * v2->x3;
v3.x3 = v1->x1 * v2->x2 - v1->x2 * v2->x1;
return(v3);
}

/* double dot_prdct(POINT *v1, POINT *v2)
{
double r;
r = v1->x1 * v2->x1 + v1->x2 * v2->x2 + v1->x3 * v2->x3;
return(r);
} */

int do_cross(POINT *Q, double *RES, double *R2, double *cos_mosaic, double *sin_mosaic)
{
  double near_rr, far_rr;
  near_rr = (*R2) - 2.0 * (*RES) * (-Q->x1 * (*cos_mosaic) + sqrt(Q->x2 * Q->x2 + Q->x3 * Q->x3) * (*sin_mosaic));
  far_rr = (*R2) - 2.0 * (*RES) * (-Q->x1 * (*cos_mosaic) - sqrt(Q->x2 * Q->x2 + Q->x3 * Q->x3) * (*sin_mosaic));
  /* Q->x1 < 0, abs(Q->x1) = -Q->x1 */
  if ((near_rr * far_rr) <= 0.0 ) {
    return (1); /* true */
  } else {
    return (0); /* false */
  }
}

int check_range(double llim, double ulim, int incld_l, int incld_u, double *value, double default_value)
{
  if (incld_l == 0) {                    /* not including lower limit */
    if (*value <= llim) {                /* including lower limit : false */
      *value = default_value;
      return (0);
    } else {
      if (incld_u == 0) {                /* not including upper limit */
	if (*value >= ulim) {            /* including upper limit : false */
	  *value = default_value;
	  return (0);
	} else {                         /* not including upper limit : true */
	  return (1);
	}
      } else {                           /* including upper limit */
	if (*value > ulim) {             /* not including upper limit : false */
	  *value = default_value;
	  return (0);
	} else {                         /* including upper limit : true */
	  return (1);
	}
      }
    }
  } else {                               /* including lower limit */
    if (*value < llim) {                 /* not including lower limit : false */
      *value = default_value;
      return (0);
    } else {
      if (incld_u == 0) {                /* not including upper limit */
	if (*value >= ulim) {            /* including upper limit : false */
	  *value = default_value;
	  return (0);
	} else {                         /* not including upper limit : true */
	  return (1);
	}
      } else {                           /* including upper limit */
	if (*value > ulim) {             /* not including upper limit : false */
	  *value = default_value;
	  return (0);
	} else {                         /* including upper limit : true */
	  return (1);
	}
      }
    }
  }
}