#define DATUM             "2010_06_15"
/****************************************************
This is the file 'mythen.c.2010_06_15', 
genuine name:    'mythen.c'.         

The aim of this software is the 2D-visualization of  the 
single exposition files made by the Mythen linear detector (1D PSD) 
e.g. at the 
"Large Four Circle X-Ray Diffractometer" and at the 
"Siemens Rotating Anode X-Ray Diffractometer" 
at the Max-Planck-Institut fuer Metallforschung (Stuttgart). 

Authors:  Agnes Szoekefalvi-Nagy, Janos Major (major@mf.mpg.de) 

The program 'mythen' is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License as published 
by the Free Software Foundation; either version 2 of the License, or (at 
your option) any later version.

This program is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
General Public License for more details.

You can obtain a copy of the GNU General Public License
from   http://www.gnu.org/copyleft/gpl.html; if not, write to 
the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
Boston, MA  02110-1301, USA.
**************************************************
   First version: mythen.c.2009_05_11
   The standard site for this file is: $HOME/c-gcc/

   To compile and link please use:    
   cp mythen.c.xxxx_xx_xx mythen.c
   cc -lm mythen.c -Wall -o ~/bin/mythen

   To make the program available for all users of the computer, please copy it 
   to the central site (usually you can do it as a root only):
   chmod 777 ~/bin/mythen
   mv ~/bin/mythen /usr/local/bin/ 

   Adjustable parameters in the lines with: */      /*P*/         /*            
**************************************************
motor conventions at the various diffractometers 
   (initial motor axes, rotation signs, positions are in [],
   v = vertical, h = horizontal, x-axis = X-ray direction):
    
#1: large four circle diffractometer, Stuttgart (fourc)
    D_MYTHEN[1]= 350.000 mm
    detector: tth [zv+]; sample: th [zv+]; 
    sample tilt: chi-90 [xh+]; sample rotation: phi [yv-: down]  

#2: Siemens rotating anode diffractometer, Stuttgart (sixc)
    D_MYTHEN[2]= 392.476 mm 
??? detector: delta [hz]; sample: th [v]; 
??? sample tilt: chi [h: X-ray]; sample rotation: phi [v: down];  
    mu = 0;   gamma = 0;    

#3: ANKA diffractometer in vertical geometry (psic)
    D_MYTHEN[3]= 450.000;
    detector: delta [zh-]; sample: eta [zh-]; 
    sample tilt: chi [xh+]; sample rotation: phi [yv+];  
    mu = 0;  nu = 0

#4: same as #3 but with standard computer

*************************************************/
#define __NO_MATH_INLINES    /* makes slower but, may be, more precise */ 

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <errno.h>

#define FORMAT "mythen <Diffractometer> <Path> <SpecFile> <ScanNo> <VerticalAxis> \
<HorizontalAxis> <Binning> <Norm>\n\n"
#define EICHUNG 0.00804 /* mythen detector: degree/channel */  /*P*/ 
#define d_MYTHEN 0.05   /* mythen detector pixel width in mm */ 
                        /* 64 mm/1280 ch = 0.05 mm/ch */
#define IMEAN 640       /* mythen detector: channel corresponding to SPEC */     
#define CHI0 90         /* basic position of the motor chi */

#define MAX_FILE 25100   /* max. number of files */
#define MAX_PIXEL 2500   /* max. number of pixels along the ppm axes */
#define MAX_CHANNEL 1400 /* max. number of detector channels*/
#define MAX_LEN 256      /* max. number of input characters */
#define EFFZERO 1e-40    /* my definition of zero */
#define NEARZERO 1e-10   /* zero for the matrix functions */
#define LIM_EXP 70.      /* limit for the function exp */
#define SCREEN_WRITE 20  /* screen write out period during interpolation */ /*P*/
#define PLOTZERO 1.e-6   /* smaller values are not plotted in gnuplot */
#define RAD 0.0174532925199433 /* 1 degree in radian */ 
#define PI 3.1415926535897932384626433832795028841971693993751

typedef struct {double x,y,z;} vector;                        /* 3D vector */
typedef struct {double xx,xy,xz,yx,yy,yz,zx,zy,zz;} matrix;   /* 3x3 tensor */

/* function prototypes: */

void message (char *text, char *text1);
void error (char *text);
void terminate (char *text);
void warten (char *text);
void read_spec_file(void);
void read_mythen_files(void);
void q_calculation (void);
void colour (FILE *file, float h, float s, float v);
void ppm (void);           /* 2D plot */
void interpolate (void);
void logarithm (void);
void grafika(void);
void write_gnuplot2D(void);
vector hkl2lab (double h, double k, double l, 
		double tth, double th, double chi, double phi);
vector angles2hkl (double tth, double th, double chi, double phi);
void distance (void);

/* function proptotypes for vector-matrix calculation: */
void orientation_matrix (void);
void printV (char *text, vector);
void printM (char *text, matrix);
vector DtimesV (double, vector);
matrix DtimesM (double D, matrix M);
double VtimesV (vector V1, vector V2);
matrix MtimesM (matrix M1, matrix M2);
vector MtimesV (matrix M, vector V);
vector VcurlV (vector V1, vector V2); 
double absV (vector V);
vector normV (vector V);
matrix MplusM (matrix M1, matrix M2);
vector VplusV (vector V1, vector V2); 
vector VminusV (vector V1, vector V2); 
matrix transpM (matrix M);
matrix Rx (double alpha);
matrix Ry (double alpha);
matrix Rz (double alpha);
matrix ROTATION (vector V, double theta);

/* global declarations: */
int diffractometer;       /* which diffractometer is used: 
                             1: large four circle diffractometer, Stuttgart 
                             2: Siemens rotating anode diffractometer, Stuttgart
                             3: ANKA diffractometer in vertical sample mode 
                                ANKA computer
                             4: ANKA diffractometer in vertical sample mode 
                                standard computer
*/
double D_MYTHEN[100];               /* possible Mythen-sample distances */  
double counts [MAX_FILE][MAX_CHANNEL];  /* experimental results */
double zeit;                            /* time of the experiment */
int MON[MAX_FILE];                          /* monitor counts */
double q_x[MAX_FILE][MAX_CHANNEL];
double q_y[MAX_FILE][MAX_CHANNEL];
double H[MAX_FILE], K[MAX_FILE], L[MAX_FILE];
double TTH[MAX_FILE], TH[MAX_FILE], CHI[MAX_FILE], PHI[MAX_FILE];
int NABS[MAX_FILE];               /* number of the absorber used */  
double eichung=EICHUNG;           /* PSD degree/channel */
double channel_m=IMEAN;           /* PSD reference channel tth=0->direct beam*/
double PPM[MAX_PIXEL][MAX_PIXEL];   /* array to be transfered to a ppm file */
double GNU[MAX_PIXEL][MAX_PIXEL];   /* array to be transferred to gnuplot */
double a_i[MAX_FILE][MAX_CHANNEL];  /* the corresponding q_x values */
double a_f[MAX_FILE][MAX_CHANNEL];  /* the corresponding alpha_j values */
char *path;                 /* path of the file names */
char *spec_file;            /* first (common) part of the file names */
char Instrument [MAX_LEN];  /* name of the diffractometer used */
char ppmfile [MAX_LEN];
char psfile [MAX_LEN];
char pngfile [MAX_LEN];
char gnuplotfile [MAX_LEN];
char gnuplot2Dfile [MAX_LEN];
int scan_Nr;
double lambda;             /* wavelength of the X-ray beam in Angstrom */
int imean=IMEAN;             /* Mythen channel number corresponding to spec */
double d_mythen=d_MYTHEN;    /* 60 mm/1280 ch = 0.05 mm/ch */ 
double D_mythen;           /* motion radius of the mean (spec) Mythen channel */
matrix O;                  /* orientation matrix */
matrix B;                  /* crystal hkl to crystal Cartesian hkl */
matrix Binv;              
double h0,k0,l0,th0,chi0,phi0,tth0;
double h1,k1,l1,th1,chi1,phi1,tth1;
float a1, a2, a3, alpha1, alpha2, alpha3;  /* real space data of the unit cell  */
float b1, b2, b3, beta1, beta2, beta3;/*reciprocal space data of the unit cell  */
vector qvn, qhn, qhor, qver;          /* plot axis */
int rlu = 1;                          /* plot units: 0 for 1/Angstrom 1 for rlu */
char Line [MAX_LEN+1];
int m_all, m_act;                 /* number of files: mythen_all, mythen_actual */
int ARGC;                         /* mode of calculation */
double absorber[100]; 
int binning = 1;                  /* number */
int CH = 1280;                    /* number of segments of the mythen detector */  
double norm = 0.0;                /* factor for the data normalization*/

/* parameters of the ppm file: */
int pixel;                        /* number of pixels along the ppm axes */      
int GEO=0;                          /* tth=th: 2 tth=th: 1 tth=const: 0 */
int LOG;                          /* 1 for log; 0 for lin */
double BG;                        /* background for the log plot */
double q_x_min;               /* range of q_x (th) */
double q_x_max;               /* range of q_x (th) */

int errorcount = 0;     /* for the function error */

/* for the functions ppm + colour: */
float colour_min = 0.79, colour_max = 0.0;                   /*P*/
int gnuplotversion;      /* 0 if >= 4.0; >0 if <=3.7 */           


/**********************    function: main    ***************/
int main(int argc, char *argv[]){
  char *kontr;
  int i;

  ARGC = argc;  
  D_MYTHEN[1]= 350.;     /* gr. 4-Kreis */                   /*P*/
  D_MYTHEN[2]= 392.476;  /* Siemens */
  D_MYTHEN[3]= 450.000;  /* ANKA vertical */
  D_MYTHEN[4]= 450.000;  /* ANKA vertical */

  for (i=0; i<100; i++) absorber[i]= 1.;
  /*
  absorber[ 0]= 1./1.000000000;
  absorber[ 1]= 1./0.551687248; 
  absorber[ 2]= 1./0.303211331;
  absorber[ 3]= 1./0.167277825;
  absorber[ 4]= 1./0.0932500738;
  absorber[ 5]= 1./0.0514448766;
  absorber[ 6]= 1./0.028274479;
  absorber[ 7]= 1./0.0155986695;
  absorber[ 8]= 1./0.00851832867;
  absorber[ 9]= 1./0.00469945331;
  absorber[10]= 1./0.00258285377;
  absorber[11]= 1./0.00142492749;
  absorber[12]= 1./0.000794334777;
  absorber[13]= 1./0.000438224368;
  absorber[14]= 1./0.000240851305;
  absorber[15]= 1./0.000132874594;
  */
  
  absorber[ 0]= 1./1.000000000;
  absorber[ 1]= 1./0.537796311; 
  absorber[ 2]= 1./0.297393608;
  absorber[ 3]= 1./0.159937185;
  absorber[ 4]= 1./0.0912335968;
  absorber[ 5]= 1./0.0490650918;
  absorber[ 6]= 1./0.0271322885;
  absorber[ 7]= 1./0.0145916447;
  absorber[ 8]= 1./0.00836996798;
  absorber[ 9]= 1./0.00450113379; 
  absorber[10]= 1./0.00248917498;
  absorber[11]= 1./0.00133866912;
  absorber[12]= 1./0.000763622284;
  absorber[13]= 1./0.000410673247;
  absorber[14]= 1./0.000227096386;
  absorber[15]= 1./0.000122131599;

  fprintf (stderr,
	   "  ____________________________________________________________________\n"

	   " |             Program version:   mythen.c." DATUM  "                 |\n"
	   " |      Max-Planck-Institut fuer Metallforschung (Stuttgart)          |\n"
	   " |                                                                    |\n"
           " | This software prepares two-dimensional plots from the data obtained|\n"
           " | by the application of the Mythen linear detector (1D PSD) at X-ray |\n"
           " | facilities. The data are normalized to 1 second or to a fraction   |\n"
           " | of the monitor counts.                                             |\n"
	   " |                                                                    |\n"
	   " | Please make sure that during the experiment the ROI of the Mythen  |\n"
           " | detector as defined by the mythen_config is centered around the    |\n"
	   " | channel imean = %3d; further if the detector motor(s) is (are) at  |\n"
	   " | zero, the detector channel imean corresponds to the incoming beam. |\n"
	   " | Please use correct value for the detector-sample distance; adjust  |\n"
	   " | if necessary.                                                      |\n"
	   " |                                                                    |\n"
	   " | The files written by the program are temporary, will be rewritten  |\n"
           " | during the next action. If you need these files, please save them. |\n"
	   " |____________________________________________________________________|\n"
	   , imean);
  switch (argc) {
    
  case 7: /* test: hkl, angles to q */
    { char ch[MAX_LEN+1];
      float h, k, l, tth, th, chi, phi;
      vector Q;
      
      GEO = 3;
      diffractometer = strtod(argv[1], &kontr);
      if(kontr == argv[1]) terminate(" wrong 1st parameter\n");
      path = argv[2];
      message("path of the data files=       ", path);
      spec_file = argv[3];
      message("             file name=", spec_file);
      scan_Nr = strtod(argv[4], &kontr);
      if(kontr == argv[4]) terminate(" wrong 4th parameter\n");
      if (scan_Nr >= MAX_FILE){
	printf("The number of files >= MAX_FILE= %i:\n", MAX_FILE);    
	terminate("please correct the C code.");}
      
      read_spec_file();
      orientation_matrix ();
      printM ("O", O);
      message("","");
      message("or0 test by hkl2lab:","");
      switch (diffractometer) {
      case 1:
        fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
		h0, k0, l0, tth0, th0, chi0+CHI0, phi0);
        break;
      case 2:
	fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
		h0, k0, l0, tth0, th0, chi0+CHI0, phi0);
        break;
      case 3:
	fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
		h0, k0, l0, tth0, th0, chi0, phi0);
        break;      
      case 4:
	fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
		h0, k0, l0, tth0, th0, chi0, phi0);
        break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }
      printV ("q = ",hkl2lab(h0,k0,l0,tth0,th0,chi0,phi0));
      message("or1 test by hkl2lab:","");
      switch (diffractometer) {
      case 1:
	fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
		h1, k1, l1, tth1, th1, chi1+CHI0, phi1);
        break;
      case 2:
	fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
		h1, k1, l1, tth1, th1, chi1+CHI0, phi1);
	break;
      case 3:
	fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
		h1, k1, l1, tth1, th1, chi1, phi1);
	break;      
      case 4:
	fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
		h1, k1, l1, tth1, th1, chi1, phi1);
	break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }
      printV ("q = ",hkl2lab(h1,k1,l1,tth1,th1,chi1,phi1));
      fprintf(stderr, "\n%s",
	      "You have started the testing mode 'hkl2lab' for a hklscan:\n");
      while (ch[0]!='q') {    
	fprintf(stderr, "%s",
		"please type 'q' for quit or 'h k l tth th chi phi':  \n");
	fgets(ch,MAX_LEN,stdin);
	sscanf (ch, "%f%f%f%f%f%f%f", &h, &k, &l, &tth, &th, &chi, &phi);
      switch (diffractometer) {
      case 1:
	chi = chi - CHI0;
	break;
      case 2:
	chi = chi - CHI0;
	break;
      case 3:
        break;      
      case 4:
        break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }
	Q.x=h; Q.y=k; Q.z=l;
	printV ("Crystal (h, k, l)             ", Q);
	Q = MtimesV (B, Q);
	printV ("Crystal Cartesian [1/Angstrom]", Q);
	Q = hkl2lab(h,k,l,tth,th,chi,phi);
	printV ("Lab [1/Angstrom]              ", Q);
      }
      return 0;
    }
    break;

  case 8: /* test: angles to hkl */
    { char ch[MAX_LEN+1];
      float ftth, fth, fchi, fphi;
      double tth,  th,  chi,  phi;
      vector Q;

      GEO = 3;
      diffractometer = strtod(argv[1], &kontr);
      if(kontr == argv[1]) terminate(" wrong 1st parameter\n");
      path = argv[2];
      message("path of the data files=       ", path);
      spec_file = argv[3];
      message("             file name=", spec_file);
      scan_Nr = strtod(argv[4], &kontr);
      if(kontr == argv[4]) terminate(" wrong 4th parameter\n");
      if (scan_Nr >= MAX_FILE){
	printf("The number of files >= MAX_FILE= %i:\n", MAX_FILE);    
	terminate("please correct the C code.");}
      
      read_spec_file();
      orientation_matrix ();
      printM ("O", O);
      message("","");
      message("or0 test by hkl2lab:","");
      switch (diffractometer) {
      case 1:
        fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
                h0, k0, l0, tth0, th0, chi0+CHI0, phi0);
        break;
      case 2:
        fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
                h0, k0, l0, tth0, th0, chi0+CHI0, phi0);
        break;
      case 3:
        fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
                h0, k0, l0, tth0, th0, chi0, phi0);
        break;      
      case 4:
        fprintf(stderr, "h0=%.2f k0=%.2f l0=%.2f tth0=%f th0=%f chi0=%f phi0=%f: \n",
                h0, k0, l0, tth0, th0, chi0, phi0);
        break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }
      printV ("q = ",hkl2lab(h0,k0,l0,tth0,th0,chi0,phi0));
      message("or1 test by hkl2lab:","");
      switch (diffractometer) {
      case 1:
        fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
                h1, k1, l1, tth1, th1, chi1+CHI0, phi1);
        break;
      case 2:
        fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
                h1, k1, l1, tth1, th1, chi1+CHI0, phi1);
        break;
      case 3:
        fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
                h1, k1, l1, tth1, th1, chi1, phi1);
        break;      
      case 4:
        fprintf(stderr, "h1=%.2f k1=%.2f l1=%.2f tth1=%f th1=%f chi1=%f phi1=%f: \n",
                h1, k1, l1, tth1, th1, chi1, phi1);
        break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }
      printV ("q = ",hkl2lab(h1,k1,l1,tth1,th1,chi1,phi1));
      message("","");
      message("or0 test by angles2hkl:","");  
      Q = angles2hkl(tth0,th0,chi0,phi0);
      fprintf(stderr, "tth0=%f th0=%f chi0=%f phi0=%f: \n"
	      "                  h0=%f k0=%f l0=%f\n",
	      tth0, th0, chi0, phi0, Q.x, Q.y, Q.z);
      message("or1 test by angles2lab:","");
      Q = angles2hkl(tth1,th1,chi1,phi1);
      fprintf(stderr, "tth1=%f th1=%f chi1=%f phi1=%f: \n"
	      "                  h1=%f k1=%f l1=%f\n",
	      tth1, th1, chi1, phi1, Q.x, Q.y, Q.z);

      fprintf(stderr, "\n%s",
	      "You have started the testing mode 'angles2hkl' for a hklscan:\n");
      while (ch[0]!='q') {    
	fprintf(stderr, "%s",
		"please type 'q' for quit or 'tth th chi phi':  \n");
	fgets(ch,MAX_LEN,stdin);
	sscanf (ch, "%f%f%f%f", &ftth, &fth, &fchi, &fphi);
      switch (diffractometer) {
      case 1:
        tth = ftth; th = fth; chi = fchi - CHI0; phi = fphi;
        break;
      case 2:
	tth = ftth; th = fth; chi = fchi - CHI0; phi = fphi;
        break;
      case 3:
	tth = ftth; th = fth; chi = fchi       ; phi = fphi;
        break;      
      case 4:
	tth = ftth; th = fth; chi = fchi       ; phi = fphi;
        break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }
      Q = angles2hkl(tth,th,chi,phi);
      fprintf(stderr, "tth=%f th=%f chi=%f phi=%f: \n"
	      "               h=%f   k=%f   l=%f\n",
	      tth, th, chi, phi, Q.x, Q.y, Q.z);
      }
      return 0;
    }
    break;

  case 6: /*horizontal axis: mythen file number 
	    vertical axis:   mythen channel number */
    GEO = 0;
    diffractometer = strtod(argv[1], &kontr);
    if(kontr == argv[1]) terminate(" wrong 1st parameter\n");
    path = argv[2];
    message("path of the data files=", path);
    spec_file = argv[3];
    message("         Spec file name=", spec_file);
    scan_Nr = strtod(argv[4], &kontr);
    if(kontr == argv[4]) terminate(" wrong 4th parameter\n");
    if (scan_Nr >= MAX_FILE){
      printf("The number of files >= MAX_FILE= %i:\n", MAX_FILE);    
      terminate("please correct the C code.");}
    binning = strtod(argv[5], &kontr);   
    if(kontr == argv[5]) terminate(" wrong 11th parameter\n");
    if (binning < 1) terminate("The binning must be at least 1.");
    break;
  
  case 13:
    GEO = 3;
    diffractometer = strtod(argv[1], &kontr);
    if(kontr == argv[1]) terminate(" wrong 1st parameter\n");
    path = argv[2];
    message("path of the data files=       ", path);
    spec_file = argv[3];
    message("             file name=", spec_file);
    scan_Nr = strtod(argv[4], &kontr);
    if(kontr == argv[4]) terminate(" wrong 4th parameter\n");
    if (scan_Nr >= MAX_FILE){
      printf("The number of files >= MAX_FILE= %i:\n", MAX_FILE);    
      terminate("please correct the C code.");}
    qver.x = strtod(argv[5], &kontr);   
    if(kontr == argv[5]) terminate(" wrong 5th parameter\n");
    qver.y = strtod(argv[6], &kontr);   
    if(kontr == argv[6]) terminate(" wrong 6th parameter\n");
    qver.z = strtod(argv[7], &kontr);   
    if(kontr == argv[7]) terminate(" wrong 7th parameter\n");
    qhor.x = strtod(argv[8], &kontr);   
    if(kontr == argv[8]) terminate(" wrong 8th parameter\n");
    qhor.y = strtod(argv[9], &kontr);   
    if(kontr == argv[9]) terminate(" wrong 9th parameter\n");
    qhor.z = strtod(argv[10], &kontr);   
    if(kontr == argv[10]) terminate(" wrong 10th parameter\n");
    binning = strtod(argv[11], &kontr);   
    if(kontr == argv[11]) terminate(" wrong 11th parameter\n");
    if (binning < 1) terminate("The binning must be at least 1.");
    norm = strtod(argv[12], &kontr);   
    if(kontr == argv[12]) terminate(" wrong 12th parameter\n");
    
    if (absV(qhor) < NEARZERO) terminate ("the vector qhor should not be so short.");
    if (absV(qver) < NEARZERO) terminate ("the vector qver should not be so short.");
    break;
    
  default:
    fprintf(stderr, "\n You have made a mistake in the input, "
	    "please try again using the format:\n\n" FORMAT
            "<Diffractometer> = 1: large four circle diffractometer, Stuttgart\n" 
            "                   2: Siemens rotating anode diffractometer, Stuttgart\n"
            "                   3: ANKA diffractometer in vertical sample mode\n"
            "                      ANKA computer \n"
            "                   4: ANKA diffractometer in vertical sample mode\n"
            "                      standard computer \n"
	    "<Path> = path of the Spec and data files (e.g., . or ../)\n"
	    "<SpecFile> = full name of the Spec file\n"
	    "<ScanNo> = Scan Number\n"
	    "<VerticalAxis>   = hkl-space coordinates defining the vertical axis \
of the plot,\n"
	    "                    e.g., 1 1 0\n"
	    "<HorizontalAxis> = hkl-space coordinates defining the horizontal axis \
of the plot,\n" 
	    "                    e.g., 0 0 1 (for an L-scan)\n"
            "<Binning> = number of channels to be combined\n" 
            "<Norm> = 0:      normalized to 1 second \n"
            "         n (>0): normalized to (monitor counts)/n \n"
	    "For the mode 'Mythen channel number vs. Mythen file number' \
please type: \n"
	    "     mythen <Diffractometer> <Path> <SpecFile> <ScanNo> <Binning>   \n"
	    "For the test mode 'h,k,l,angles -> q' please type: \n"
	    "     mythen <Diffractometer> <Path> <SpecFile> <ScanNo> x x   \n"
	    "For the test mode 'angles -> h,k,l' please type: \n"
	    "     mythen <Diffractometer> <Path> <SpecFile> <ScanNo> x x x  \n"
	    "The C source code of this program is available from:\n"
	    "http://www.mf.mpg.de/en/abteilungen/dosch/software/software_en.shtml\n");
    fprintf(stderr, "The correct format for the usage of this program is:\n" FORMAT);
    terminate ("  input error  ");
    break;}
  
  sprintf(gnuplot2Dfile, "%s%s%s%i%s", spec_file,"_", "scan",scan_Nr, ".2D");
  
  /* the gnuplot version will be asked: */
  {FILE *pipe;
    pipe = popen("gnuplot", "w");
    fprintf(pipe, "set terminal png \n");
    gnuplotversion =  pclose(pipe);
    /*if (gnuplotversion == 0) message("A new gnuplot version is present.","\n");
    else message
      ("An old gnuplot version is present: no png files will be available.",
      "\n");*/}
  /* ready: gnuplot version */  
  switch (diffractometer) {
  case 1: strcpy (Instrument, "<large 4-c Stgt>"); break;
  case 2: strcpy (Instrument, "<Siemens Stgt>"); break;
  case 3: strcpy (Instrument, "<ANKA vertical>"); break;
  case 4: strcpy (Instrument, "<ANKA vertical>"); break;
  default:  terminate ("illegal diffractometer number.\n" FORMAT);}
  read_spec_file();
  orientation_matrix ();
  read_mythen_files ();
  q_calculation ();  
  write_gnuplot2D();
  grafika();
  error ("terminal"); 
  return 0;}
/**********************    function: main    ***************/

/**********************    function: hkl2lab    ***************/
vector hkl2lab (double h, double k, double l, 
		double tth, double th, double chi, double phi) {
  /* transforms the crystal (h,k,l) vector to the laboratory coordinates
     x: initial beam direction
     y: horizontal  
     z: vertical downwards                                */
  vector q;
  q.x = h; q.y = k; q.z = l;
  q = MtimesV (B, q);
  q = MtimesV (O ,q);
  q = MtimesV (Ry(phi*RAD), q);
  q = MtimesV (Rx(chi*RAD), q);
  q = MtimesV (Rz( th*RAD), q);
  return q;}
/**********************    function: hkl2lab    ***************/

/**********************  function: angles2hkl   ***************/
vector angles2hkl (double tth, double th, double chi, double phi) {
  /* transforms the angular positions to the crystal (h,k,l) vector */
  vector q;
  double k00 = 2*PI/lambda;
  q.x = k00*(-1.+cos(tth*RAD)); q.y = k00*sin(tth*RAD); q.z = 0.;
  q = MtimesV (Rz( -th*RAD), q);
  q = MtimesV (Rx(-chi*RAD), q);
  q = MtimesV (Ry(-phi*RAD), q);
  q = MtimesV (transpM(O),   q);
  q = MtimesV (Binv,         q);
  return q;}
/**********************  function: angles2hkl   ***************/

/**********************    function: read_spec_file    ***************/
void read_spec_file () {
  int i;
  char *kontr;
  char line [MAX_LEN+1],*token;
  char dummy [MAX_LEN], d [100][MAX_LEN]; 
  char motor [MAX_LEN] ;
  float start, stop, step;  
  char motor1 [MAX_LEN],motor2[MAX_LEN];
  float start1, start2, stop1, stop2;
  float fH, fK, fL;
  float fh1, fk1, fl1, fth1, fchi1, fphi1, ftth1;
  float fh0, fk0, fl0, fth0, fchi0, fphi0, ftth0, flambda;
  float fa1, fa2, fa3, falpha1, falpha2, falpha3;
  float fb1, fb2, fb3, fbeta1, fbeta2, fbeta3;
  float ftth, fth, fchi, fphi, fzeit;
  int idummy, n_attn=0, n_mon=0; 
  FILE *file; 
  char p_spec_file [MAX_LEN];
  int scan_Nr_found = 0;

  sprintf (p_spec_file, "%s/%s", path, spec_file);


  file=fopen(p_spec_file,"r");
  if (file==NULL)  terminate("the desired spec file does not exist.\n" FORMAT);
  while (fgets (Line, MAX_LEN, file) != NULL) {
    if ((Line[0]=='#')&&(Line[1]=='S')){
      sscanf (Line, "%s%i%s ", 
	      dummy,&idummy,dummy);
      if (idummy==scan_Nr) {
	scan_Nr_found = 1;
	printf("scan: %s ", Line);
	printf("\nscanNr.=: %i     scan=:  %s     ",scan_Nr,dummy);
	if (strcmp(dummy, "ascan")); else {GEO = 0.;
	  sscanf (Line, "%s%s%s%s%f%f%f%f", 
		  dummy, dummy, dummy,motor,&start,&stop,&step,&fzeit);
	  m_all=step+1;}
	if (strcmp(dummy, "a2scan")); else {GEO = 0;
	  sscanf (Line, "%s%s%s%s%f%f%s%f%f%f%f", 
		  dummy, dummy, dummy,motor1,&start1,&stop1,
		  motor2,&start2,&stop2,&step,&fzeit);
	  m_all=step+1;
	  
	  printf("m_all=%i     step=%f             motor1=%s      start1=%f "
		 "  stop1=%f       motor2=%s        start2=%f stop2=%f \n",
		 m_all,step,motor1,start1,stop1,motor2,start2,stop2) ;}
	if (strcmp(dummy, "hklscan")); else {	
	  sscanf (Line, "%s%s%s%s%s%s%s%s%s%f%f", 
		  dummy, dummy, dummy,dummy, dummy, dummy,dummy, 
		  dummy, dummy,&step,&fzeit);
	  m_all=step+1;}
	break;     }   } }

  if (scan_Nr_found == 0) terminate("wrong scan number. \n" FORMAT);

  while (fgets (line, MAX_LEN, file) != NULL) {
    if ((line[0]=='#')&&(line[1]=='G')&&(line[2]=='1')){
      sscanf (line, 
	      "%s%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%f%s%s%f%f%f%f%s%s%f%s", 
	      dummy,
	      &fa1, &fa2, &fa3, &falpha1, &falpha2, &falpha3,
	      &fb1, &fb2, &fb3, &fbeta1,  &fbeta2,  &fbeta3,
	      &fh0, &fk0, &fl0, &fh1, &fk1, &fl1, 
	      &ftth0, &fth0, &fchi0, &fphi0, dummy, dummy,
	      &ftth1, &fth1, &fchi1, &fphi1, dummy, dummy, &flambda, dummy);
      switch (diffractometer) {
      case 1:
	tth0=ftth0; th0=fth0;         chi0=fchi0-CHI0;         phi0=-fphi0;
	tth1=ftth1; th1=fth1;         chi1=fchi1-CHI0;         phi1=-fphi1;
	break;
      case 2:
	tth0=ftth0; th0=fth0;         chi0=fchi0-CHI0;         phi0=-fphi0;
	tth1=ftth1; th1=fth1;         chi1=fchi1-CHI0;         phi1=-fphi1;
	break;
      case 3:
	tth0=-ftth0; th0=-fth0; chi0=fchi0; phi0=fphi0;
	tth1=-ftth1; th1=-fth1; chi1=fchi1; phi1=fphi1;
	break;      
      case 4:
	tth0=-ftth0; th0=-fth0; chi0=fchi0; phi0=fphi0;
	tth1=-ftth1; th1=-fth1; chi1=fchi1; phi1=fphi1;
	break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }

      a1=fa1; a2=fa2; a3=fa3; alpha1=falpha1; alpha2=falpha2; alpha3=falpha3;
      b1=fb1; b2=fb2; b3=fb3;  beta1=fbeta1;   beta2=fbeta2;   beta3=fbeta3;
      lambda = flambda;
      zeit = fzeit;
      
      /*printf("scan=%i \nh0=%f k0=%f l0=%f tth0=%f th0=%f chi0=%f phi0=%f"
	"\nh1=%f k1=%f l1=%f tth1=%f th1=%f chi1=%f phi1=%f lambda=%f zeit=%f\n",
	scan_Nr,h0,k0,l0,tth0,th0,chi0,phi0,
	h1,k1,l1,tth1,th1,chi1,phi1,lambda,zeit); */
      break;     }   }
  
  while (fgets (line, MAX_LEN, file) != NULL) {
    if ((line[0]=='#')&&(line[1]=='L')){
      sscanf (line,"%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s",
	      d[ 0],d[ 1],d[ 2],d[ 3],d[ 4],d[ 5],d[ 6],d[ 7],d[ 8],d[ 9],
              d[10],d[11],d[12],d[13],d[14],d[15],d[16],d[17],d[18],d[19],
              d[20],d[21],d[22],d[23],d[24],d[25],d[26],d[27],d[28],d[29],
              d[30]);
      /*fprintf (stderr, "d = %s \n",d);*/
      for (i=1;i<=30;i++)
	if ((d[i][0]=='a')&&(d[i][1]=='t')&&(d[i][2]=='t')&&(d[i][3]=='n')) 
	  {n_attn=i; break;}
      for (i=1;i<=30;i++)
	if ((d[i][0]=='M')&&(d[i][1]=='o')&&(d[i][2]=='n')&&(d[i][3]=='i')) 
	  {n_mon =i; break;}
      break;     }}
  
  if (n_attn>0) fprintf(stderr,
			"Attenuators may have been applied during the experiment:\n"
			"            %2ith column of the spec file. \n", n_attn);
  else  fprintf(stderr,
		"Attenuators have not been applied during the experiment. \n");
  if (n_mon>0) fprintf(stderr,
		       "The monitor counts are in the "
		       "%2ith column of the spec file. \n", n_mon);
  
  m_act=0;

  while (fgets (line, MAX_LEN, file) != NULL) {
    if ((line[0]=='#')||(strlen(line)==1))  break;
    sscanf (line,  
	    "%f%f%f%f%f%f%f%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s", 
	    &fH, &fK, &fL, &ftth, &fth, &fchi, &fphi, d[ 8],d[ 9],
	    d[10],d[11],d[12],d[13],d[14],d[15],d[16],d[17],d[18],d[19],
	    d[20],d[21],d[22],d[23],d[24],d[25],d[26],d[27],d[28],d[29],
	    d[30]);
    
    m_act=m_act+1;  
    H[m_act]    = fH;
    K[m_act]    = fK;
    L[m_act]    = fL;
    if (n_mon>0) MON[m_act] = strtod(d[n_mon],&kontr); else MON[m_act] = 0;
    if (MON[m_act] == 0) norm = 0;

      switch (diffractometer) {
      case 1:
	TTH[m_act]  = ftth; 
	TH[m_act]  = fth;
	CHI[m_act] = fchi-CHI0;    
	PHI[m_act] = -fphi;
	break;
      case 2:
	TTH[m_act]  = ftth; 
	TH[m_act]  = fth;
	CHI[m_act] = fchi-CHI0;    
	PHI[m_act] = -fphi;
	break;
      case 3: 
	TTH[m_act]  = -ftth;
	TH[m_act]  = -fth;
	CHI[m_act] = fchi;    
	PHI[m_act] = fphi;
	break;      
      case 4: 
	TTH[m_act]  = -ftth;
	TH[m_act]  = -fth;
	CHI[m_act] = fchi;    
	PHI[m_act] = fphi;
	break;      
      default:  terminate ("illegal diffractometer number.\n" FORMAT); }

      if (n_attn>0) NABS[m_act]= strtod(d[n_attn],&kontr); else NABS[m_act]= 0;
    /*fprintf(stderr, "absorber Nr.= %i factor = %f\n", 
      NABS[m_act], absorber[NABS[m_act]]);*/ } 

  /* B matrix transforms from the crystal h,k,l to crystal Cartesian q.x,q.y,q.z:
     W. R. Busing, H. A. Levy, Appl. Cryst. (1967) 22, 457 */

  B.xx = b1;
  B.xy = b2*cos(beta3*RAD);
  B.xz = b3*cos(beta2*RAD);
  B.yx = 0;
  B.yy = b2*sin(beta3*RAD);
  B.yz =-b3*sin(beta2*RAD)*cos(alpha1*RAD);
  B.zx = 0;
  B.zy = 0;
  B.zz = 2.*PI/a3;

  Binv.xx = 1./b1;
  Binv.xy = -cos(beta3*RAD)/sin(beta3*RAD)/b1;
  Binv.xz = a3/(2.*PI*b1*b2*sin(beta3*RAD))*
		(-b2*b3*cos(alpha1*RAD)*cos(beta3*RAD)*sin(beta2*RAD)
		-b2*b3*cos(beta2*RAD)*sin(beta3*RAD));
  Binv.yx = 0.;
  Binv.yy = 1./(b2*sin(beta3*RAD));
  Binv.yz = a3*b3*cos(alpha1*RAD)*sin(beta2*RAD)/(2.*PI*b2*sin(alpha1*RAD));
  Binv.zx = 0.;
  Binv.zy = 0.;
  Binv.zz = a3/(2.*PI);

  printM("B",B); 

  h0=fh0; k0=fk0; l0=fl0; 
  h1=fh1; k1=fk1; l1=fl1;

  /* the newline character will be removed (because of gnuplot): */
  token = strtok(Line, "\n");
  strcpy(Line, token);

  /* printf ("kanal = %i   %i  \n", kanal, intens); */     
  fclose(file);
  printf("m_all=%i  m_act=%i \n", m_all, m_act) ;
  return;}
/**********************    function: read_spec_file    ***************/

/**********************    function: read_mythen_files    ***************/
void read_mythen_files () {
  /* experimental counts will be read in:  */
  int i, ii, k, kanal, intens;
  char filename [MAX_LEN];
  char line [MAX_LEN+1];
  double sum;
  FILE *file; 

  for (ii=1; ii<=m_act; ii++){
    sprintf(filename, "%s/%s%s%s%i%s%i%s%i%s", path, spec_file,"_", "scan",
	    scan_Nr,"_",ii,"of", m_all, ".myt");
    
    file=fopen(filename,"r");
    if (file==NULL)  { 
      fprintf(stderr, "The desired file %s does not exist.\n",filename);
      terminate(FORMAT);}
    while (fgets (line, MAX_LEN, file) != NULL) {
      sscanf (line, "%i%i", &kanal, &intens);

      if (norm == 0) 
	counts[ii][kanal]=intens * absorber [NABS[ii]] / zeit;
      else 
	counts[ii][kanal]=intens * absorber [NABS[ii]] / (double) MON[ii] * norm;}
    fclose(file);}
  
  CH = CH / binning;
  for (ii=1; ii<=m_act; ii++){
    for (i=0; i<CH; i++) {
      sum = 0;
      for (k=i*binning; k<(i+1)*binning; k++)
	sum = sum + counts[ii][k];
      counts[ii][i]=sum;}}
  eichung  =  eichung * binning;
  d_mythen = d_mythen * binning;
  imean    = imean / binning;
  if (norm == 0)
    fprintf (stderr, "The mythen counts are normalized to 1 second.\n");
  else
    fprintf (stderr, 
	     "The mythen counts are normalized to (monitor counts)/%6.2g.\n",
	     norm);
  return;}
/**********************    function: read_mythen_files    ***************/

/**********************    function: q_calculation    ***************/
void q_calculation(){
/* momentum vector  will be calculated to every experimental counts: */
  int ii,jj;
  vector q, q0;
  double k00;
  double delta;
  
  D_mythen = D_MYTHEN[diffractometer];

  switch (GEO){

   case 0:                                  
    /* horizontal: number of the mythen file  
       vertical: mythen-channel number    */
    for (ii=1; ii<=m_act; ii++)
      for (jj=0; jj<=CH-1; jj++){                           
	q_x[ii][jj] = ii;
	q_y[ii][jj] = jj;}          /*   or	q_y[ii][jj] = CH-1-jj     */   
    break;

   case 1:                                    /* tth = th */
    /* reflectometry-sample horizon: the same channel */
    for (ii=1; ii<=m_act; ii++)
      for (jj=0; jj<=CH-1; jj++){                           
	q_x[ii][jj] = q_x_min + 
	  (q_x_max - q_x_min)/(m_act-1.) * (ii-1.);
	q_y[ii][jj] = (jj-channel_m)*eichung;}     
    break;
    
   case 2:                                    /* tth = 2 * th */
    /* reflectometry-specular direction: the same channel */
    for (ii=1; ii<=m_act; ii++)
      for (jj=0; jj<=CH-1; jj++){                           
	q_x[ii][jj] = q_x_min + 
	  (q_x_max - q_x_min)/(m_act-1.) * (ii-1.);
	q_y[ii][jj] = q_x[ii][jj] + (jj-channel_m)*eichung;}     
    break;

  case 3:                                          /* diffractometry */
    k00=2*PI/lambda;
    
    for (ii=1; ii<=m_act; ii++){
      qhn = hkl2lab(qhor.x,qhor.y,qhor.z,TTH[ii],TH[ii],CHI[ii],PHI[ii]);
      if (rlu == 0) qhn = normV(qhn);
      else qhn = DtimesV(1./absV(qhn)/absV(qhn),qhn);  
      qvn = hkl2lab(qver.x,qver.y,qver.z,TTH[ii],TH[ii],CHI[ii],PHI[ii]);
      if (rlu == 0)  qvn = normV(qvn); 
      else qvn = DtimesV(1./absV(qvn)/absV(qvn),qvn);  
      q0  = hkl2lab(H[ii], K[ii], L[ii], TTH[ii],TH[ii],CHI[ii],PHI[ii]);
      
      for (jj=0; jj<=CH-1; jj++){    /* actual Mythen geometry :*/
	delta = atan(d_mythen*(jj-imean)/D_mythen); 
	/*fprintf(stderr,"delta = %f \n",delta/RAD);*/
	q.x = k00*(cos(delta)-1.)+q0.x*cos(delta)-q0.y*sin(delta);
	q.y = q0.y*cos(delta)+(q0.x+k00)*sin(delta);
        q.z = q0.z;
	
	if ((ii==1)||(ii==m_act/2)||(ii==m_act))
	  if ((jj==0)||(jj==imean)||(jj==CH-1)) {
	    printf("file No=%4i  channel No=%4i q.x=%10.4g q.y=%10.4g q.z=%10.4g ",
		   ii, jj, q.x, q.y, q.z);
	    if (jj==imean) printf("*\n"); else printf("\n");}
	
	q_x[ii][jj]= VtimesV(qhn,q);     
	q_y[ii][jj]= VtimesV(qvn,q);} }
    break;
    
  default: terminate ("geometry parameter (GEO): illegal value ");
    break;}
  
  return;}
/**********************    function: q_calculation    ***************/

/***********************  write_gnuplot2D******************/
void  write_gnuplot2D(void){
  FILE *file; 
  int i,j;
  time_t curtime;
  
 /* makes 2D gnuplot file */
    /* against compiler warning message: */
  file = fopen(gnuplot2Dfile, "w"); fclose( file);  
  curtime = time (NULL);
  file = fopen (gnuplot2Dfile, "w") ;
  if ((ARGC == 13)&&(diffractometer == 3))
    fprintf (file, 
	     "#  \n"
	     "# This file was produced by the software 'mythen' version %s \n"
	     "# %s  \n" 
	     "# %s  \n"
	     "# %s  \n"
	     "# D_mythen = %6.3f mm \n"
	     "# Instrument: %s \n"
	     "# %s"
	     "# Binning: %i \n"
	     "# Norm: %6.2g (0: for 1 s; >0: for (monitor count)/value)\n"
	     "# Components of the momentum transfer:   \n"
	     "#       first  column: along the hkl-direction (%f,%f,%f)    \n"    
	     "#       second column: along the hkl-direction (%f,%f,%f)    \n"  
	     "# Normalized counts:  \n" 
	     "#       third  column    \n"  
	     "#  \n", 
	     DATUM, path, spec_file ,Line, D_MYTHEN[diffractometer], Instrument,  
	     asctime (localtime (&curtime)), binning, norm, 
	     qhor.x,qhor.y,qhor.z,qver.x,qver.y,qver.z);

  if ((ARGC == 13)&&(diffractometer != 3))
    fprintf (file, 
	     "#  \n"
	     "# This file was produced by the software 'mythen' version %s \n"
	     "# %s  \n" 
	     "# %s  \n"
	     "# %s  \n"
	     "# D_mythen = %6.3f mm \n"
	     "# Instrument: %s \n"
	     "# %s"
	     "# Binning: %i \n"
	     "# Norm: %6.2g (0: for 1 s; >0: for (monitor count)/value)\n"
	     "# Components of the momentum transfer:   \n"
	     "#       first  column: along the hkl-direction (%f,%f,%f)    \n"    
	     "#       second column: along the hkl-direction (%f,%f,%f)    \n"  
	     "# Normalized counts:  \n" 
	     "#       third  column    \n"  
	     "#  \n", 
	     DATUM, path, spec_file ,Line, D_MYTHEN[diffractometer], Instrument, 
	     asctime (localtime (&curtime)), binning, norm,
	     qhor.x,qhor.y,qhor.z,qver.x,qver.y,qver.z);
  
  if (ARGC == 6)
    fprintf (file, 
	     "#  \n"
	     "# This file was produced by the software 'mythen' version %s \n"
	     "# %s  \n" 
	     "# %s  \n"
	     "# %s  \n"
	     "# D_mythen = %6.3f mm \n"
	     "# Instrument: %s \n"
	     "# %s"
	     "# Binning: %i \n"
	     	     "# First  column: Mythen file number    \n"    
	     "# Second column: Mythen channel number    \n"  
	     "# Third  column: detected counts/second    \n"  
	     "#  \n", 
	     DATUM, path, spec_file ,Line, D_MYTHEN[diffractometer], Instrument, 
	     asctime (localtime (&curtime)), binning);
  
  for (i=1; i<=m_act; i++){ 
    for (j=0; j<=CH-1; j++) 
      fprintf (file, " %f \t  %f \t  %f  \n", q_x[i][j], q_y[i][j], counts[i][j]);
    fprintf (file, " \n");} 
  fclose(file);
  printf ("\nThe file \"%s\" has been written or rewritten. \n\n",gnuplot2Dfile);}
/***********************  write_gnuplot2D******************/

/********************  function: grafika *******************/
/* GNUPLOT HOWTO:               t '<titel>' w l 
   (in newer gnuplot versions:) lw <linewidth> lt <colour> 
   (in older gnuplot versions:) lw <linewidth> <colour> 
   w e == with errorbars */
/* specular case: equation: q_y = m * q_x
   distance: sqrt(1+m^2)*(q_y-m*q_x)
   all other:     equation: q_x/q_x0 + q_y/q_y0 = 1
   distance:
   (q_x0^2+q_y0^2)*(q_x0*q_y0-q_x0*q_y-q_y0*q_x)
*/

void grafika(void){
  char ch[MAX_LEN+1];
  FILE *pipe;  
  int logscale = 1;  
  int plotscaler = 0;  
  float CL, CH;
  
  fprintf (stderr, "Sample-detector distance = %6.3f mm \n\n", 
	   D_MYTHEN[diffractometer]);
  
  pipe = popen("gnuplot -geometry 600x600+10+10", "w");
  if (diffractometer != 3) fprintf (pipe, "set term wxt noraise \n"); 
  fprintf (pipe, "set colorbox user origin .87,.15 size .04,.7 \n");
 start:
  if (logscale == 1) {fprintf (pipe, "set logscale z \n");
    fprintf (pipe, "set logscale cb \n");}
  else {fprintf (pipe, "set nologscale z \n");
    fprintf (pipe, "set nologscale cb \n");}
  fprintf (pipe, "set view 0,0 \n");
  fprintf (pipe, "set grid \n");
  if (GEO == 3) { 
    if (rlu==0) {
      fprintf (pipe, "set title \"%s        scan No: %4i\
          D_mythen=%.3f mm  \\n %s  %s \\n y-axis: (%.0f  %.0f  %.0f)	 [Angstrom^-1]\" \n", 
	       spec_file, scan_Nr, D_mythen, Instrument, Line, 
	       qver.x, qver.y, qver.z);
      fprintf (pipe, "set xlabel 'x-axis: (%.0f  %.0f  %.0f) [Angstrom^-1]' \n", 
	       qhor.x, qhor.y, qhor.z ); } 
    else {
      fprintf (pipe, "set title \"%s        scan No: %4i\
          D_mythen=%.3f mm   \\n %s  %s \\n y-axis: (%.0f  %.0f  %.0f)	[rlu]\" \n",
	       spec_file, scan_Nr, D_mythen, Instrument, Line, 
	       qver.x, qver.y, qver.z);
      fprintf (pipe,
	       "set xlabel 'x-axis: (%.0f  %.0f  %.0f) [rlu]' \n",
	       qhor.x, qhor.y, qhor.z ); } }
  else {
    fprintf (pipe, "set title \"%s       scan No: %4i\
	     D_mythen=%.3f mm   \\n  %s  %s \\n y-axis: Mythen channel number \" \n", 
	     spec_file, scan_Nr, D_mythen, Instrument, Line);
    fprintf (pipe, 
	     "set xlabel 'x-axis: Mythen file number' \n"); }
  fprintf (pipe, "set ticslevel 0 \n" );
  fflush(pipe);
  
  fprintf (pipe,    
	   "splot [:][:][:] '%s' t '' with pm3d \n", gnuplot2Dfile);
  fflush(pipe); 
  fprintf(stderr, 
	  "'i' = linear scale for counts   'o' = logarithmic scale for counts\n"
	  "'r' = reciprocal-lattice units  'a' = Angstrom^-1 unit \n"
	  "'g' = refresh the graphics      "
          "'d' = adjust the detector-sample distance \n"
	  "'c' = rescale the colour scale         \n"
	  "'p' = (re)write graphic files (ps/eps/png/pdf) \n"     
	  "'q' = quit                      Further with RETURN: "); 
  fgets(ch,MAX_LEN,stdin);
  if (ch[0]=='i') {logscale=0; goto start;} 
  else
    if (ch[0]=='o') {logscale=1; goto start;} 
    else
      if (ch[0]=='r') 
	{rlu=1; q_calculation(); write_gnuplot2D(); goto start;} 
      else
	if (ch[0]=='a') 
	  {rlu=0; q_calculation(); write_gnuplot2D(); goto start;} 
	else
	  if (ch[0]=='g') {fclose (pipe);  
	    pipe = popen("gnuplot -geometry 600x600+10+10", "w");
	    if (diffractometer != 3) fprintf (pipe, "set term wxt noraise \n");     
	    fprintf (pipe, "set colorbox user origin .87,.15 size .04,.7 \n");
	    goto start;} 
	  else
	    if (ch[0]=='p') {        
	      plotscaler = plotscaler + 1;

	      fprintf(pipe, "set terminal postscript color solid lw 2\n");     
	      fprintf(pipe, "set output '%s'\n", "mythen.ps");
	      fprintf(pipe, "replot \n");
	      
	      fprintf(pipe, "set terminal postscript eps color solid lw 2\n"); 
    	      fprintf(pipe, "set output '%s' \n", "mythen.eps");
	      fprintf(pipe, "replot \n");

	      fprintf(pipe, "set terminal png size 1280, 960 \n");      
              fprintf(pipe, "set output '%s' \n", "mythen.png"); 
	      fprintf(pipe, "replot \n");

	      if (diffractometer != 3) {
		fprintf(pipe, "set terminal jpeg size 1280, 960 \n");    
		fprintf(pipe, "set output '%s' \n", "mythen.jpg");    
		fprintf(pipe, "replot \n"); }                           
	      
	      if (diffractometer == 3) {
		fprintf(pipe, "set terminal pdf size 6, 6 \n");    
		fprintf(pipe, "set output '%s' \n", "mythen.pdf"); 
		fprintf(pipe, "replot \n"); }                       
	      
	      goto start;}
	    else 
	      if (ch[0]=='d') {distance (); goto start;}
	      else 
		if (ch[0]=='c') {  
		  fprintf(stderr, 
			  "Please type the new lower and higher limits"
			  " for the colour scale:");
		  fgets(ch, MAX_LEN, stdin);
		  sscanf (ch, "%f %f ", &CL, &CH); 
		  fprintf (pipe,"set cbrange [%f:%f] \n", CL, CH);
		  goto start;}
		else
		  if ((ch[0]=='q')||(ch[0]=='Q')) ;
		  else goto start; 
  fclose (pipe);}
/***********************  function: grafika ************************/

/***********************  function: distance ************************/
void distance (void) {/* Changes the predefined distance between the 
                         sample centre and the Mythen detector (channel 
                         imean) == motion radius of the mean (imean: 
                         spec) Mythen channel. */
  char ch [MAX_LEN+1];
  float fd;
  
  fprintf(stderr,
	  "The distance between the sample centre and the Mythen detector\n"
	  "is defined as: %6.3f mm.\n"
	  "Please type: RETURN if you agree with the present value or\n"
	  "             the desired distance in mm and RETURN ---> ",
	  D_MYTHEN[diffractometer]);
  fgets(ch,MAX_LEN,stdin);
  sscanf (ch, "%f", &fd);
  if (fd > NEARZERO) D_MYTHEN[diffractometer] = fd;
  
  fprintf (stderr, "The new sample-detector distance = %6.3f mm \n\n", 
	   D_MYTHEN[diffractometer]); 
  q_calculation(); write_gnuplot2D(); 
  return;}
/***********************  function: distance ************************/

/********************** function: orientation_matrix ***************/
void orientation_matrix () {
  /* Calulates the orientation matrix O of the sample. */
  
  vector s0, s1; /*momentum transfer calculated from the angles */
  vector t0, t1; /*momentum transfer calculated from h, k, l*/
  vector d0, d1; /* difference vectors */
  vector ki;
  vector r; /* rotation axis of O*/
  double theta, theta0, theta1; /* rotation angle of O */
  vector dummyV, dummyV0, dummyV1, dummyVs, dummyVt;
  double dummy, k00;

  k00=2*PI/lambda;
  ki.x = k00; ki.y = 0; ki.z = 0;

  dummyV.x =  -1+cos(tth0*RAD);
  dummyV.y =     sin(tth0*RAD);
  dummyV.z =   0;
  printV("q0",dummyV);
  dummyV = MtimesV (Rz(- th0*RAD),dummyV);
  dummyV = MtimesV (Rx(-chi0*RAD),dummyV);
  dummyV = MtimesV (Ry(-phi0*RAD),dummyV);
  s0 = normV(dummyV);
  dummyV.x = h0;
  dummyV.y = k0;
  dummyV.z = l0;
  t0 = normV(MtimesV(B,dummyV));
  d0=VminusV(s0,t0);

  dummyV.x =  -1+cos(tth1*RAD);
  dummyV.y =     sin(tth1*RAD);
  dummyV.z =   0          ;
  printV("q1",dummyV);
  dummyV = MtimesV (Rz(- th1*RAD),dummyV);
  dummyV = MtimesV (Rx(-chi1*RAD),dummyV);
  dummyV = MtimesV (Ry(-phi1*RAD),dummyV);
  s1 = normV(dummyV);
  dummyV.x = h1;
  dummyV.y = k1;
  dummyV.z = l1;
  t1 = normV(MtimesV(B,dummyV));
  d1=VminusV(s1,t1);

  /* calculation of the rotation axis r:  */
  if (absV(VcurlV(d0,d1)) > NEARZERO) {/* normal case */
    fprintf (stderr, "Calculation of the orientation matrix: normal case\n");
    r = normV(VcurlV(d0,d1));}
  else { /* d0 and d1 are parallel NOT YET TESTED */
    fprintf (stderr, "Calculation of the orientation matrix: special case\n");
    dummyV0 = VplusV(DtimesV(1./absV(d0),s0),DtimesV(1./absV(d0),t0));
    dummyV1 = VplusV(DtimesV(1./absV(d1),s1),DtimesV(1./absV(d1),t1));
    dummyV  = VminusV(dummyV0, dummyV1);
    r = normV(dummyV);};
  fprintf (stderr, "The rotation axis of the orientation matrix:\n");
  printV("r",r);
 /* calculation of the rotation angle theta: */

  dummyVs = DtimesV(VtimesV(s0,r),r);
  dummyVs = VminusV(s0,dummyVs);
  dummyVt = DtimesV(VtimesV(t0,r),r);
  dummyVt = VminusV(t0,dummyVt);
  dummy = VtimesV(dummyVs, dummyVt)/absV(dummyVs)/absV(dummyVt);
  theta0 = acos(dummy);            /* in radian unit */

  if (VtimesV(dummyVs,r)>1.e-6) printf("nem meroleges a: s \n");
  if (VtimesV(dummyVt,r)>1.e-6) printf("nem meroleges a: t \n");

  dummyVs = DtimesV(VtimesV(s1,r),r);
  dummyVs = VminusV(s1,dummyVs);
  dummyVt = DtimesV(VtimesV(t1,r),r);
  dummyVt = VminusV(t1,dummyVt);
  dummy = VtimesV(dummyVs, dummyVt)/absV(dummyVs)/absV(dummyVt);
  theta1 = acos(dummy);           /* in radian unit */

  if (VtimesV(dummyVs,r)>1.e-6) printf("nem meroleges b: s \n");
  if (VtimesV(dummyVt,r)>1.e-6) printf("nem meroleges b: t \n");

  theta = (theta0 + theta1)/2.;  /* these numbers should be the same */
  
  O = ROTATION (r, theta);
  fprintf (stderr, "The rotation angle of the orientation matrix:\n");
  printf("    theta0: %g  theta1: %g     mean: %g [degrees]  \n", 
	 theta0/RAD, theta1/RAD, theta/RAD);
  return;}
/********************** function: orientation_matrix ***************/

/********************** function: printV ***************/
void printV (char *text, vector V) {
  /* prints out the elements of the vector V */
  fprintf(stderr, "%20s: .x =%10.4g  .y =%10.4g   .z =%10.4g   \n", 
	  text, V.x, V.y, V.z);
  return;}
/********************** function: printV ***************/

/********************** function: printM ***************/
void printM (char *text, matrix M) {
  /*  prints out the elements of the vector M */
  fprintf(stderr, "%20s:    .xx=%10.4g  .xy=%10.4g   .xz=%10.4g  \n "
	          "%19s     .yx=%10.4g  .yy=%10.4g   .yz=%10.4g  \n "
	          "%19s     .zx=%10.4g  .zy=%10.4g   .zz=%10.4g  \n ",
	  text, M.xx, M.xy, M.xz, "", M.yx, M.yy, M.yz, "", M.zx, M.zy, M.zz);
  return;}
/********************** function: printM ***************/

/********************** function: DtimesV ***************/
vector DtimesV (double D, vector V) {
  /*  multiplies the vector V by the scalar D */
  V.x = V.x * D;
  V.y = V.y * D;
  V.z = V.z * D;
  return (V);}
/********************** function: DtimesV ***************/

/********************** function: DtimesM ***************/
matrix DtimesM (double D, matrix M) {
  /*  multiplies the matrix M by the scalar D */
  M.xx = M.xx * D;
  M.xy = M.xy * D;
  M.xz = M.xz * D;
  M.yx = M.yx * D;
  M.yy = M.yy * D;
  M.yz = M.yz * D;
  M.zx = M.zx * D;
  M.zy = M.zy * D;
  M.zz = M.zz * D;
  return (M);}
/********************** function: DtimesM ***************/

/********************** function: VtimesV ***************/
double VtimesV (vector V1, vector V2) {
  /*  skalar product of the vectors V1 and V2 */
  return (V1.x*V2.x+V1.y*V2.y+V1.z*V2.z);}
/********************** function: VtimesV ***************/

/********************** function: MtimesM ***************/
matrix MtimesM (matrix M1, matrix M2) {
  /*  multiplies the matrices M1 and M2  */
  matrix M;
  M.xx = M1.xx*M2.xx+M1.xy*M2.yx+M1.xz*M2.zx;
  M.xy = M1.xx*M2.xy+M1.xy*M2.yy+M1.xz*M2.zy;
  M.xz = M1.xx*M2.xz+M1.xy*M2.yz+M1.xz*M2.zz;
  M.yx = M1.yx*M2.xx+M1.yy*M2.yx+M1.yz*M2.zx;
  M.yy = M1.yx*M2.xy+M1.yy*M2.yy+M1.yz*M2.zy;
  M.yz = M1.yx*M2.xz+M1.yy*M2.yz+M1.yz*M2.zz;
  M.zx = M1.zx*M2.xx+M1.zy*M2.yx+M1.zz*M2.zx;
  M.zy = M1.zx*M2.xy+M1.zy*M2.yy+M1.zz*M2.zy;
  M.zz = M1.zx*M2.xz+M1.zy*M2.yz+M1.zz*M2.zz;
  return (M);}
/********************** function: MtimesM ***************/

/********************** function: MtimesV ***************/
vector MtimesV (matrix M, vector V) {
  /*  multiplies the matrix M and the vector V  */
  vector VR;
  VR.x=M.xx*V.x+M.xy*V.y+M.xz*V.z;
  VR.y=M.yx*V.x+M.yy*V.y+M.yz*V.z;
  VR.z=M.zx*V.x+M.zy*V.y+M.zz*V.z;
  return (VR);}
/********************** function: MtimesV ***************/

/********************** function: VcurlV ***************/
vector VcurlV (vector V1, vector V2) {
  /*  vector product of the vectors V1 and V2 */
  vector VR;
  VR.x=V1.y*V2.z-V1.z*V2.y;
  VR.y=V1.z*V2.x-V1.x*V2.z;
  VR.z=V1.x*V2.y-V1.y*V2.x;
  return (VR);}
/********************** function: VcurlV ***************/

/********************** function: absV ***************/
double absV (vector V) {
  /*  absolute value (magnitude) of the vector V */
  return (sqrt(V.x*V.x+V.y*V.y+V.z*V.z));}
/********************** function: absV ***************/

/********************** function: normV ***************/
vector normV (vector V) {
  /*  makes a unity vector of the vector V */
  double D;
  D=absV(V);
  if (D<EFFZERO) 
    terminate ("normV: the vector is a null vector - cannot proceed");
  return (DtimesV(1./D,V));}
/********************** function: normV ***************/

/********************** function: MplusM ***************/
matrix MplusM (matrix M1, matrix M2) {
  /*  sums of the matrices M1 and M2  */
  matrix M;
  M.xx = M1.xx+M2.xx;
  M.xy = M1.xy+M2.xy;
  M.xz = M1.xz+M2.xz;
  M.yx = M1.yx+M2.yx;
  M.yy = M1.yy+M2.yy;
  M.yz = M1.yz+M2.yz;
  M.zx = M1.zx+M2.zx;
  M.zy = M1.zy+M2.zy;
  M.zz = M1.zz+M2.zz;
  return (M);}
/********************** function: MplusM ***************/

/********************** function: VplusV ***************/
vector VplusV (vector V1, vector V2) {
  /*  sums of the vectors V1 and V2  */
  vector V;
  V.x = V1.x+V2.x;
  V.y = V1.y+V2.y;
  V.z = V1.z+V2.z;
  return (V);}
/********************** function: VplusV ***************/

/********************** function: VminusV ***************/
vector VminusV (vector V1, vector V2) {
  /*  difference of the vectors V1 and V2  */
  vector V;
  V.x = V1.x-V2.x;
  V.y = V1.y-V2.y;
  V.z = V1.z-V2.z;
  return (V);}
/********************** function: VminusV ***************/

/********************** function: transpM ***************/
matrix transpM (matrix M) {
  /*  transposes of the matrix M  */
  /*  for rotation matrices this is the inverse */
  matrix MR;
  MR.xx = M.xx;
  MR.xy = M.yx;
  MR.xz = M.zx;
  MR.yx = M.xy;
  MR.yy = M.yy;
  MR.yz = M.zy;
  MR.zx = M.xz;
  MR.zy = M.yz;
  MR.zz = M.zz;
  return (MR);}
/********************** function: transpM ***************/

/********************** function: Rx ***************/
matrix Rx (double alpha) {
  /* rotational matrix around the x axis*/
  /* alpha argument in radian */
  matrix M;
  double c,s;
  c = cos(alpha);
  s = sin(alpha);
  M.xx = 1;
  M.xy = 0;
  M.xz = 0;
  M.yx = 0;
  M.yy = c;
  M.yz = -s;
  M.zx = 0;
  M.zy = s;
  M.zz = c;
  return M;}
/********************** function: Rx ***************/

/********************** function: Ry ***************/
matrix Ry (double alpha) {
  /* rotational matrix around the y axis*/
  /* alpha argument in radian */
  matrix M;
  double c,s;
  c =  cos(alpha);
  s =  sin(alpha);   		    
  M.xx = c;
  M.xy = 0;
  M.xz = s;
  M.yx = 0;
  M.yy = 1;
  M.yz = 0;
  M.zx = -s;
  M.zy = 0;
  M.zz = c;
  return M;}
/********************** function: Ry ***************/

/********************** function: Rz ***************/
matrix Rz (double alpha) {
  /* rotational matrix around the z axis*/
  /* alpha argument in radian */
  matrix M;
  double c,s;
  c = cos(alpha);
  s = sin(alpha);
  M.xx = c;
  M.xy = -s;
  M.xz = 0;
  M.yx = s;
  M.yy = c;
  M.yz = 0;
  M.zx = 0;
  M.zy = 0;
  M.zz = 1;
  return M;}
/********************** function: Rz ***************/

/********************** function: ROTATION ***************/
matrix ROTATION (vector V, double theta) {
  /* V is the axis of the rotation
     theta is the angle of rotation in radian
     This function calculates the corresponding rotation matrix */
  matrix M, P, I, Q;
  P.xx = V.x*V.x;
  P.xy = V.x*V.y;
  P.xz = V.x*V.z;
  P.yx = V.y*V.x;
  P.yy = V.y*V.y;
  P.yz = V.y*V.z;
  P.zx = V.z*V.x;
  P.zy = V.z*V.y;
  P.zz = V.z*V.z;
  I.xx = 1.;
  I.xy = 0.;
  I.xz = 0.;
  I.yx = 0.;
  I.yy = 1.;
  I.yz = 0.;
  I.zx = 0.;
  I.zy = 0.;
  I.zz = 1.;
  Q.xx =   0.;
  Q.xy = -V.z;
  Q.xz =  V.y;
  Q.yx =  V.z;
  Q.yy =   0.;
  Q.yz = -V.x;
  Q.zx = -V.y;
  Q.zy =  V.x;
  Q.zz =   0.;
  M = MplusM(DtimesM(cos(theta),I),DtimesM(sin(theta),Q));
  M = MplusM(DtimesM(1.-cos(theta),P),M);
  return M;}
/********************** function: ROTATION ***************/



/********************    function: terminate    ***************/
void terminate (char *text){                   /* Text bei Fehler */
message("terminate:", text);
exit (1);}
/********************    function: terminate    ***************/

/********************      function: message    ***************/
void message (char *text1, char *text2){        /*Text beim Suchen*/
fprintf(stderr, "%s %s \n", text1, text2);}
/********************      function: message    ***************/

/********************      function: error    ***************/
void error (char *ref){                /* reports on runtime errors 
                                          Usage:       error ("label");
                          #include <fenv.h> and #include <errno.h> and 
                          int errorcount = 0;  (global) are necessary.*/
#ifdef _FENV_H   
  int raised;  
#endif /* fenv.h*/
  char text[MAX_LEN];  

  strcpy(text, "errno error message (");
  strcat(text, ref);  strcat(text, ")");
  perror(text);

#ifdef _FENV_H 
  raised = fetestexcept 
    (FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
  /*(FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);*/
  if (raised) {errorcount = errorcount + 1; 
  sprintf(text, "%-3i", errorcount);
  strcat(text, "fp error message (");
  strcat(text, ref);  strcat(text, "):");
  /*if (raised & FE_INEXACT)   strcat(text, " INEXACT ");*/  
  if (raised & FE_DIVBYZERO) strcat(text, " DIVBYZERO ");
  if (raised & FE_UNDERFLOW) strcat(text, " UNDERFLOW ");
  if (raised & FE_OVERFLOW)  strcat(text, " OVERFLOW ");
  if (raised & FE_INVALID)   strcat(text, " INVALID ");
  /* putchar('\a'); fflush(stdout); */
  message (text, " ");    
  whistle();}
  /* if (raised) terminate("floating point error"); */
  feclearexcept (FE_ALL_EXCEPT);    /**/  
#endif /* fenv.h*/
  errno = 0; /**/                  }
/******************     function: error    ***************/

/********************    function: warten    ******************/
void warten (char *text){
  char ch[MAX_LEN];  
  fprintf(stderr, "Further with RETURN:  %s", text);
  fgets(ch, MAX_LEN-1, stdin);}
/********************    function: warten    ******************/

/*********************  function: whistle  *******************/
void whistle (void){putchar('\a'); fflush(stdout);}
/********************** function: whistle  *******************/