/* crysconv.f -- translated by f2c (version 20000817).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
  doublereal s[16] /* was [4][4] */ , sinv[16] /* was [4][4] */ ;
  integer iok;
} osolem_;

#define osolem_1 osolem_

/* 	================ */

/* dec$	  Ident	    'V01A' */
/* ------------------------------------------------------------------ */
/* Updates: */
/*  V01A  7-May-1996 DM.	  Put error output to IOLUN, use IMPLICIT NONE and */
/* 			    get the code indented so that it is readable! */
/*  made F77 compliant. Mark Koennecke July 1996 */
/* ------------------------------------------------------------------ */
/*   Routines to deal with the reciprocical lattice PB */
/* ------------------------------------------------------------------ */
/*   Entry points in this file: */

/*       SETRLP : CALCULATION OF S AND INVS , ORIENTATION MATRIX */
/*       RL2SPV : TRANSFO FROM RECIP LAT TO SCAT PLANE */
/*       SP2RLV : TRANSFO FROM SCAT PLANE TO RECIP LAT */
/* 	INVS   : INVERT MATRIX S, GENERATED BY SETRLP. */
/*       ERRESO : DEAL ITH ERROR MESSAGES FOR ALL MODULES */

/*        SUBROUTINE SETRLP(SAM,IER) */
/*        SUBROUTINE RL2SPV(QHKL,QT,QM,QS,IER) */
/*        SUBROUTINE SP2RLV(QHKL,QT,QM,QS,IER) */
/* 	 SUBROUTINE INVS(S,SINV,IER) */
/*        SUBROUTINE ERRESO(MODULE,IER) */
/* ------------------------------------------------------------------ */
/* Subroutine */ int setrlp_(doublereal * sam, integer * ier)
{
  /* System generated locals */
  doublereal d__1;

  /* Builtin functions */
  double cos(doublereal), sin(doublereal), sqrt(doublereal),
      atan(doublereal);

  /* Local variables */
  static doublereal alfa[3], cosa[3], cosb[3];
  static integer imod;
  static doublereal sina[3], sinb[3], aspv[6] /* was [3][2] */ ;
  extern /* Subroutine */ int invs_(doublereal *, doublereal *, integer *);
  static doublereal a[3], b[3], c__[3], bb[9] /* was [3][3] */ , cc;
  static integer id, ie, jd, je, jf, kg, lf, lh, md, me, ne;
  static doublereal zp, vv[9] /* was [3][3] */ ;
  extern /* Subroutine */ int erreso_(integer *, integer *);
  static doublereal rlb[6] /* was [3][2] */ ;

/* 	============================ */

/*  SETRLP:    Computation of matrix S which transforms (QH,QK,QL) to */
/*       vector (Q1,Q2) in scattering plane (defined by vectors A1,A2) */
/*       and SINV matrix for the inverse transformation */

/*       INPUT SAM  SAMPLE CHARACTERISTICS */
/*       SAM(1)=AS       LATTICE PARAMETERS */
/*       SAM(2)=BS       ------------------ */
/*       SAM(3)=CS       ------------------ */
/*       SAM(4)=AA       LATTICE ANGLES */
/*       SAM(5)=BB       -------------- */
/*       SAM(6)=CC       -------------- */
/*       SAM(7)=AX       VECTOR A IN SCATTERING PLANE */
/*       SAM(8)=AY       ---------------------------- */
/*       SAM(9)=AZ       ---------------------------- */
/*       SAM(10)=BX      VECTOR B IN SCATTERING PLANE */
/*       SAM(11)=BY      ---------------------------- */
/*       SAM(12)=BZ      ---------------------------- */
/*       OUTPUT  IER ERROR RETURN TO BE TREATED BY ERRESO */
/*       IER=1   ERROR ON LATTICE PARAMETERS */
/*       IER=2   ERROR ON LATTICE ANGLES */
/*       IER=3   ERROR ON VECTORS A1, A2 */
/* ------------------------------------------------------------------ */

/* ------------------------------------------------------------------ */
/*     			  Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* 	DO NOT EXPORT THE FOLLOWING COMMON ! */
/*       IT IS JUST FOR PERMANENT STORAGE USE */

/* ----------------------------------------------------------------------- */

/* ----------------------------------------------------------------------- */
/*       SOME TESTS AND INIT OF CALCUALTION */

  /* Parameter adjustments */
  --sam;

  /* Function Body */
  *ier = 0;
  imod = 1;
  zp = 6.2831853071795862;
  osolem_1.iok = 0;
  for (id = 1; id <= 3; ++id) {
    a[id - 1] = sam[id];
    alfa[id - 1] = sam[id + 3];
    aspv[id - 1] = sam[id + 6];
    aspv[id + 2] = sam[id + 9];
/* L10: */
  }

  for (id = 1; id <= 3; ++id) {
    *ier = 1;
    if ((d__1 = a[id - 1], abs(d__1)) <= 1e-8) {
      goto L999;
    }
    *ier = 0;
/* L20: */
  }
  for (id = 1; id <= 3; ++id) {
    a[id - 1] /= zp;
    alfa[id - 1] /= 57.29577951308232087679815481410517;
    cosa[id - 1] = cos(alfa[id - 1]);
    sina[id - 1] = sin(alfa[id - 1]);
/* L30: */
  }
  cc = cosa[0] * cosa[0] + cosa[1] * cosa[1] + cosa[2] * cosa[2];
  cc = cosa[0] * 2. * cosa[1] * cosa[2] + 1. - cc;
  *ier = 2;
  if (cc <= .1) {
    goto L999;
  }
  *ier = 0;
  cc = sqrt(cc);
  je = 2;
  kg = 3;
  for (id = 1; id <= 3; ++id) {
    b[id - 1] = sina[id - 1] / (a[id - 1] * cc);
    cosb[id - 1] =
        (cosa[je - 1] * cosa[kg - 1] -
         cosa[id - 1]) / (sina[je - 1] * sina[kg - 1]);
    sinb[id - 1] = sqrt(1. - cosb[id - 1] * cosb[id - 1]);
    rlb[id + 2] = (d__1 = atan(sinb[id - 1] / cosb[id - 1]), abs(d__1)) *
        57.29577951308232087679815481410517;
    je = kg;
    kg = id;
/* L40: */
  }
  bb[0] = b[0];
  bb[1] = 0.;
  bb[2] = 0.;
  bb[3] = b[1] * cosb[2];
  bb[4] = b[1] * sinb[2];
  bb[5] = 0.;
  bb[6] = b[2] * cosb[1];
  bb[7] = -b[2] * sinb[1] * cosa[0];
  bb[8] = 1. / a[2];

  for (id = 1; id <= 3; ++id) {
    rlb[id - 1] = 0.;
    for (je = 1; je <= 3; ++je) {
/* Computing 2nd power */
      d__1 = bb[je + id * 3 - 4];
      rlb[id - 1] += d__1 * d__1;
/* L60: */
    }
    *ier = 1;
    if ((d__1 = rlb[id - 1], abs(d__1)) <= 1e-8) {
      goto L999;
    }
    *ier = 0;
    rlb[id - 1] = sqrt(rlb[id - 1]);
/* L50: */
  }
/* ----------------------------------------------------------------------- */
/* GENERATION OF S ORIENTATION MATRIX REC. LATTICE TO SCATTERING PLANE */

  for (kg = 1; kg <= 2; ++kg) {
    for (ie = 1; ie <= 3; ++ie) {
      vv[kg + ie * 3 - 4] = 0.;
      for (jf = 1; jf <= 3; ++jf) {
        vv[kg + ie * 3 - 4] += bb[ie + jf * 3 - 4] * aspv[jf + kg * 3 - 4];
/* L90: */
      }
/* L80: */
    }
/* L70: */
  }
  for (md = 3; md >= 2; --md) {
    for (ne = 1; ne <= 3; ++ne) {
      id = md % 3 + 1;
      je = (md + 1) % 3 + 1;
      kg = ne % 3 + 1;
      lh = (ne + 1) % 3 + 1;
      vv[md + ne * 3 - 4] = vv[id + kg * 3 - 4] * vv[je + lh * 3 - 4] -
          vv[id + lh * 3 - 4] * vv[je + kg * 3 - 4];
/* L110: */
    }
/* L100: */
  }

  for (id = 1; id <= 3; ++id) {
    c__[id - 1] = 0.;
    for (je = 1; je <= 3; ++je) {
/* Computing 2nd power */
      d__1 = vv[id + je * 3 - 4];
      c__[id - 1] += d__1 * d__1;
/* L130: */
    }
    *ier = 3;
    if ((d__1 = c__[id - 1], abs(d__1)) <= 1e-6) {
      goto L999;
    }
    *ier = 0;
    c__[id - 1] = sqrt(c__[id - 1]);
/* L120: */
  }

  for (id = 1; id <= 3; ++id) {
    for (je = 1; je <= 3; ++je) {
      vv[je + id * 3 - 4] /= c__[je - 1];
/* L160: */
    }
/* L150: */
  }
  for (kg = 1; kg <= 3; ++kg) {
    for (me = 1; me <= 3; ++me) {
      osolem_1.s[kg + (me << 2) - 5] = 0.;
      for (lf = 1; lf <= 3; ++lf) {
        osolem_1.s[kg + (me << 2) - 5] += vv[kg + lf * 3 - 4] * bb[lf
                                                                   +
                                                                   me * 3 -
                                                                   4];
/* L190: */
      }
/* L180: */
    }
/* L170: */
  }
  osolem_1.s[15] = 1.;
  for (jd = 1; jd <= 3; ++jd) {
    osolem_1.s[(jd << 2) - 1] = 0.;
    osolem_1.s[jd + 11] = 0.;
/* L200: */
  }
/* ----------------------------------------------------------------------- */
/*       INVERT TRANSFORMATION MATRIX S AND PU RESULT IN SINV */

  *ier = 3;
  invs_(osolem_1.s, osolem_1.sinv, ier);
  *ier = 0;
  if (*ier != 0) {
    goto L999;
  }
  osolem_1.iok = 123;
/* --------------------------------------------------------------------------- */
/*       SORTIE */

L999:
  if (*ier != 0) {
    erreso_(&imod, ier);
  }
  return 0;
}                               /* setrlp_ */

/* =========================================================================== */
/* Subroutine */ int rl2spv_(doublereal * qhkl, doublereal * qt,
                             doublereal * qm,
                             doublereal * qs, integer * ier)
{
  /* System generated locals */
  doublereal d__1;

  /* Builtin functions */
  double sqrt(doublereal);

  /* Local variables */
  static integer id, je;

/* 	========================================= */

/* ------------------------------------------------------------------ */
/*       INPUT QHKL		      QHKL -> QT */
/* 	A Q-VECTOR TO BE TRANSFORM FROM RECIP LATTICE TO SCATTERING PLANE */
/*       CHECK THAT Q-VECTOR IS IN THE PLANE */

/*       INPUT Q-VECTOR  QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE */

/*       OUTPUT Q-VECTOR QT(3) Q-VECTOR IN SCATTERING PLANE */
/*       OUTPUT QM AND QS  QMODULUS AND ITS SQUARE ( TO BE  VERIFIED ) */
/*       OUTPUT ERROR IER */
/*       IER=1   MATRIX S NOT OK */
/*       IER=2   Q NOT IN SCATTERING PLANE */
/*       IER=3   Q MODULUS TOO SMALL */
/* ------------------------------------------------------------------ */

/* ------------------------------------------------------------------ */
/* 					    Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* 	DO NOT EXPORT THE FOLLWING COOMON ! */
/*       IT IS JUST FOR PERMANENT STORAGE USE */

/* --------------------------------------------------------------------------- */
/* --------------------------------------------------------------------------- */
/* 	INIT AND TEST IF TRANSFO MATRICES ARE OK */
  /* Parameter adjustments */
  --qt;
  --qhkl;

  /* Function Body */
  *ier = 1;
  if (osolem_1.iok != 123) {
    goto L999;
  }
  *ier = 0;
/* ----------------------------------------------------------------------- */

  for (id = 1; id <= 3; ++id) {
    qt[id] = 0.;
    for (je = 1; je <= 3; ++je) {
      qt[id] += qhkl[je] * osolem_1.s[id + (je << 2) - 5];
/* L20: */
    }
/* L10: */
  }
  *ier = 2;
  if (abs(qt[3]) > 1e-4) {
    goto L999;
  }
  *ier = 0;
  *qs = 0.;
  for (id = 1; id <= 3; ++id) {
/* Computing 2nd power */
    d__1 = qt[id];
    *qs += d__1 * d__1;
/* L30: */
  }
  if (*qs < 1e-8) {
    *ier = 3;
  } else {
    *qm = sqrt(*qs);
  }
/* --------------------------------------------------------------------------- */

L999:
  return 0;
}                               /* rl2spv_ */

/* =========================================================================== */
/* Subroutine */ int sp2rlv_(doublereal * qhkl, doublereal * qt,
                             doublereal * qm,
                             doublereal * qs, integer * ier)
{
  /* System generated locals */
  doublereal d__1;

  /* Builtin functions */
  double sqrt(doublereal);

  /* Local variables */
  static integer id, je;

/* 	========================================= */

/* ------------------------------------------------------------------ */
/*       INPUT QT	   QHKL <- QT */
/* A Q-VECTOR TO BE TRANSFORM FROM SCATTERING PLANE TO RECIP LATTICE */
/*       CHECK THAT Q, D & G VECTORS ARE IN THE SCATTERING PLANE */

/*       INPUT Q-VECTOR  QT(3) Q-VECTOR IN SCATTERING PLANE */

/*       OUTPUT Q-VECTOR QHKL(3) Q-VECTOR IN RECIPROCICAL LATTICVE */
/*       OUTPUT QM AND QS  QMODULUS AND ITS SQUARE ( TO BE  VERIFIED ) */
/*       OUTPUT ERROR IER */
/*       IER=1   MATRIX S NOT OK */
/*       IER=2   Q NOT IN SCATTERING PLANE */
/*       IER=3   Q MODULUS TOO SMALL */
/* ------------------------------------------------------------------ */

/* ------------------------------------------------------------------ */
/* 					    Define the dummy arguments */
/* ------------------------------------------------------------------ */
/* 	DO NOT EXPORT THE FOLLWING COOMON ! */
/*       IT IS JUST FOR PERMANENT STORAGE USE */

/* --------------------------------------------------------------------------- */
/* --------------------------------------------------------------------------- */
/* 	INIT AND TEST IF TRANSFO MATRICES ARE OK */
  /* Parameter adjustments */
  --qt;
  --qhkl;

  /* Function Body */
  *ier = 1;
  if (osolem_1.iok != 123) {
    goto L999;
  }
  *ier = 2;
  if (abs(qt[3]) > 1e-4) {
    goto L999;
  }
  *ier = 0;
/* ----------------------------------------------------------------------- */
  *qs = 0.;
  for (id = 1; id <= 3; ++id) {
/* Computing 2nd power */
    d__1 = qt[id];
    *qs += d__1 * d__1;
/* L10: */
  }
  if (*qs < 1e-8) {
    *ier = 3;
  } else {
    *qm = sqrt(*qs);
  }
/* ----------------------------------------------------------------------- */

  for (id = 1; id <= 3; ++id) {
    qhkl[id] = 0.;
    for (je = 1; je <= 3; ++je) {
      qhkl[id] += osolem_1.sinv[id + (je << 2) - 5] * qt[je];
/* L30: */
    }
/* L20: */
  }
/* --------------------------------------------------------------------------- */

L999:
  return 0;
}                               /* sp2rlv_ */

/* ========================================================================== */
/* Subroutine */ int invs_(doublereal * s, doublereal * sinv,
                           integer * ier)
{
  /* Initialized data */

  static integer m[3] = { 2, 3, 1 };
  static integer n[3] = { 3, 1, 2 };

  static integer id, je, mi, mj, ni, nj;
  static doublereal det;

/* 	============================== */

/* ------------------------------------------------------------------ */
/*       ROUTINE TO INVERT MATRIX S, GENERATED BY SETRLP, WHICH TRANSFORMS */
/*       (QH,QK,QL) TO (Q1,Q2) IN THE SCATTERING PLANE */
/*       INPUT MATRIX DOUBLE PRECISION S(4,4) */
/*       OUTPUT MATRIX DOUBLE PRECISION SINV(4,4) */
/*       OUTPUT ERROR IER */
/*       IER=1   DETERMINANT OF MATRIX S TOO SMALL */
/* ------------------------------------------------------------------ */

/* ------------------------------------------------------------------ */
/* 					    Define the dummy arguments */
/* ------------------------------------------------------------------ */

  /* Parameter adjustments */
  sinv -= 5;
  s -= 5;

  /* Function Body */
/* ------------------------------------------------------------------ */
  *ier = 0;
  for (id = 1; id <= 4; ++id) {
    for (je = 1; je <= 4; ++je) {
      sinv[id + (je << 2)] = 0.;
/* L20: */
    }
/* L10: */
  }
  det = 0.;
  for (id = 1; id <= 3; ++id) {
    for (je = 1; je <= 3; ++je) {
      mi = m[id - 1];
      mj = m[je - 1];
      ni = n[id - 1];
      nj = n[je - 1];
      sinv[je + (id << 2)] =
          s[mi + (mj << 2)] * s[ni + (nj << 2)] - s[ni +
                                                    (mj << 2)] * s[mi +
                                                                   (nj <<
                                                                    2)];
/* L40: */
    }
    det += s[id + 4] * sinv[(id << 2) + 1];
/* L30: */
  }
  if (abs(det) < 1e-6) {
    *ier = 1;
  } else {
    for (id = 1; id <= 3; ++id) {
      for (je = 1; je <= 3; ++je) {
        sinv[id + (je << 2)] /= det;
/* L70: */
      }
/* L60: */
    }
  }
  sinv[20] = 1.;
  return 0;
}                               /* invs_ */

/* ========================================================================= */
/* Subroutine */ int erreso_(integer * module, integer * ier)
{

  /* System generated locals */
  integer i__1;

  /* Local variables */
  static integer lmod, lier;

/* 	============================= */

/* ------------------------------------------------------------------ */
/*       SUBROUTINE TO TREAT ERRORS FROM RESOLUTION CALCULATIONS */
/*     	  MODULE = 1 -> SETRLP */
/*     	  MODULE = 2 -> RL2SPV */
/*     	  MODULE = 3 -> EX_CASE */
/* ------------------------------------------------------------------ */


/* 	INCLUDE 'MAD_DEF:IOLSDDEF.INC' */
/* ------------------------------------------------------------------ */
/* 					    Define the dummy arguments */
/* ------------------------------------------------------------------ */


/* --------------------------------------------------------------------------- */
/* Computing MIN */
  i__1 = max(*ier, 1);
  lier = min(i__1, 4);
/* Computing MIN */
  i__1 = max(*module, 1);
  lmod = min(i__1, 3);
/* 	WRITE(6,501) MESER(LIER,LMOD) */
/*  501	FORMAT(A) */
  return 0;
}                               /* erreso_ */