#define DATUM                  "2009_06_02"
/****************************************************
This is the file 'offspecular.c.2009_06_02', 
genuine name:    'offspecular.c'.         

The aim of this software is to combine the single exposition files made
by the Mythen linear detector (1D PSD) at the 
"Large Four Circle X-Ray Diffractometer" at the 
Max-Planck-Institut fuer Metallforschung (Stuttgart) in reflectivity
mode. In the first part a graphics (ppm) file is constructed, in the
second part different cuts along straight lines in the alpha_i-alpha_f
plane in the form of graphics files (ps, png).  

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

The program 'offspecular' 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: offspecular.c.2009_01_27
   The standard site for this file is: $HOME/c-gcc/

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

   Adjustable parameters in the lines with: */      /*P*/         /*            
                                  
*************************************************/
 
#define __NO_MATH_INLINES    /* makes slower but, may be, more precise */ 

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

#define EICHUNG 0.00804 /* mythen detector: degree/channel */  /*P*/
#define MAX_FILE 5100    /* max. number of files */
#define MAX_PIXEL 2500   /* max. number a 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 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 */

/* function prototypes: */

void message (char *text, char *text1);
void error (char *text);
void terminate (char *text);
void wait (char *text);
void read_files(void);
void angles (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);
double rocking(double alpha_i0,double alpha_f0,double alpha_i,double alpha_f);
double specular(double m, double alpha_i, double  alpha_f);  
double horizontal(double alpha_f0, double alpha_f);
double vertical(double alpha_i0, double alpha_i);

/* global declarations: */
int    counts [MAX_FILE][MAX_CHANNEL];  /* experimental results */
double alpha_i[MAX_FILE][MAX_CHANNEL];
double alpha_f[MAX_FILE][MAX_CHANNEL];
double eichung=EICHUNG;                 /* PSD degree/channel */
double channel_m;                 /* 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 transfered to gnuplot */
double a_i[MAX_FILE][MAX_CHANNEL];  /* the corresponding alpha_i values */
double a_f[MAX_FILE][MAX_CHANNEL];  /* the corresponding alpha_j values */
int number;            /* number of files */
char *name;            /* first (common) part of the file names */
char ppmfile [MAX_LEN];
char psfile [MAX_LEN];
char pngfile [MAX_LEN];
char gnuplotfile [MAX_LEN];
char gnuplot2Dfile [MAX_LEN];

/* parameters of the ppm file: */
int pixel;                        /* number a pixels along the ppm axes */      
int GEO;                          /* 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 alpha_i_min;               /* range of alpha_i (th) */
double alpha_i_max;               /* range of alpha_i (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;
  message ("\nProgram version:   offspecular.c." DATUM, "" );
  if (argc != 10) {message("\n You have made a mistake in the input,",
                          "please try again using the format:\n");
    message (" ", 
 "offspecular <file><number><geometry><log><pixel><a_i_min><a_i_max><ch_m><BG>\n");
    message ("<file> =","first (common) part of the file names");
    message ("<number> =","number of files");
    message ("<geometry> =","tth=2th: 2  tth=th: 1  tth=const: 0");
    message ("<log> =","log10: 1  lin: 0");
    message ("<pixel> =",
	     "number of pixels in the ppm file (horizontal and vertical)");
    message ("<a_i_min> =","th_min");
    message ("<a_i_max> =","th_max");
    message ("<ch_m> =",
	     "the mean mythen channel, e.g. specular or sample horizon");
    message ("<BG> =","background for the log plot");
    message ("","");
    terminate ("  input error  ");}
  
  name = argv[1];
  message("filename=", name);
  
  number = strtod(argv[2], &kontr);
  if(kontr == argv[2]) terminate(" wrong second parameter\n");
  
  if (number >= MAX_FILE){
    printf("The number of files >= MAX_FILE= %i:\n", MAX_FILE);
    terminate("please correct the C code.");}

  GEO = strtod(argv[3], &kontr);
  if(kontr == argv[3]) terminate(" wrong third parameter\n");

  LOG = strtod(argv[4], &kontr);
  if(kontr == argv[4]) terminate(" wrong fourth parameter\n");

  pixel = strtod(argv[5], &kontr);
  if(kontr == argv[5]) terminate(" wrong fifth parameter\n");

  if (pixel >= MAX_PIXEL){
    printf("The number of pixel >= MAX_PIXEL= %i:\n", MAX_PIXEL);
    terminate("please correct the C code or use a smaller pixel number.");}

  alpha_i_min = strtod(argv[6], &kontr);
  if(kontr == argv[6]) terminate(" wrong sixth parameter\n");

  alpha_i_max = strtod(argv[7], &kontr);
  if(kontr == argv[7]) terminate(" wrong seventh parameter\n");

  channel_m = strtod(argv[8], &kontr);
  if(kontr == argv[8]) terminate(" wrong eighth parameter\n");

  BG = strtod(argv[9], &kontr);
  if(kontr == argv[9]) terminate(" wrong nineth parameter\n");

  sprintf(ppmfile, "%s%s", name, ".ppm");
  sprintf(gnuplotfile, "%s%s", name, ".gnu");
  sprintf(gnuplot2Dfile, "%s%s", name, ".2D");
  sprintf(psfile, "%s%s", name, ".ps");
  sprintf(pngfile, "%s%s", name, ".png");
  
  /* 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 */  
 
  read_files ();
  message ("\n\nProgram version:   offspecular.c." DATUM, "\n" );
  printf("\nFirst part: the 2D ppm file will be constructed:\n"
	 "\n     For reflectivity scans:"     
	 "\n         horizontal axis: grazing initial angle "
	 "\n         vertical   axis: grazing final   angle "
	 "\n         both axes have the same angular range. \n");
  angles (); 
  interpolate();
  if (LOG == 1) logarithm ();
  ppm ();  
  printf("Second part: the 1D cuts will be constructed:\n\n");
  grafika ();
  error ("terminal"); 
  return 0;}
/**********************    function: main    ***************/

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

  for (ii=1; ii<=number; ii++){
    sprintf(filename, "%s%i%s%i%s", name, ii, "of", number, ".myt");
    printf("%s \n", filename);
    
    file=fopen(filename,"r");
    if (file==NULL)  terminate(" The desired file does not exist.  ");
    while (fgets (line, MAX_LEN, file) != NULL) {
      sscanf (line, "%i%i", &kanal, &intens);
      counts[ii][kanal]=intens;  
      /* printf ("kanal = %i   %i  \n", kanal, intens); */ }
    fclose(file);}
  
  return;}
/**********************    function: read_files    ***************/


/**********************    function: angles    ***************/
void angles(){
/* physical angles will be calculated to every experimental counts: */
  int ii,jj;

  switch (GEO){

  case 2:                                    /* tth = 2 * th */
    /* specular direction: the same channel */
    for (ii=1; ii<=number; ii++)
      for (jj=0; jj<=1279; jj++){                           
	alpha_i[ii][jj] = alpha_i_min + 
	  (alpha_i_max - alpha_i_min)/(number-1.) * (ii-1.);
	alpha_f[ii][jj] = alpha_i[ii][jj] + (jj-channel_m)*eichung;}     
    break;
    
  case 1:                                    /* tth = th */
    /* sample horizon: the same channel */
    for (ii=1; ii<=number; ii++)
      for (jj=0; jj<=1279; jj++){                           
	alpha_i[ii][jj] = alpha_i_min + 
	  (alpha_i_max - alpha_i_min)/(number-1.) * (ii-1.);
	alpha_f[ii][jj] = (jj-channel_m)*eichung;}     
    break;
    
  case 0:                                   /* tth = constant */
    /* rocking geometry */
    for (ii=1; ii<=number; ii++)
      for (jj=0; jj<=1279; jj++){                           
	alpha_i[ii][jj] = alpha_i_min + 
	  (alpha_i_max - alpha_i_min)/(number-1.) * (ii-1.);
	alpha_f[ii][jj] = - alpha_i[ii][jj] + (jj-channel_m)*eichung;}     
    break;
    
  default: terminate ("geometry parameter (GEO): illegal value ");
    break;}
  
  return;}
/**********************    function: angles    ***************/


/**********************    function: interpolate    ***************/
void interpolate(){
  /* method: IDW -- Inverse Distance Weighting 
     the matrix PPM[][] will be calculated. 
delta = (max-min)/(pixel-1)
raster (PPM): alpha_i= min+(i-1)*delta
              alpha_f= min+(j-1)*delta
              i = 1,..., pixel
              j = 1,..., pixel 
experiment:   alpha_i=alpha_i[ii][jj]
              alpha_f=alpha_f[ii][jj]
              ii = 1,..., number
              jj = 0,..., 1279 
*/
  
  int i,j, ii, jj;
  double sigma, sigma2,d2, di, df, w, sum, norm, plot_min, plot_max, delta, arg; 

  /* determination of the smallest and largest angles: */
  plot_min = alpha_i[1][0];
  plot_max = alpha_i[1][0];
  for (ii=1; ii<=number; ii++)
    for (jj=0; jj<=1279; jj++) {
      if (plot_min > alpha_i[ii][jj]) plot_min = alpha_i[ii][jj];
      if (plot_min > alpha_f[ii][jj]) plot_min = alpha_f[ii][jj];
      if (plot_max < alpha_i[ii][jj]) plot_max = alpha_i[ii][jj];
      if (plot_max < alpha_f[ii][jj]) plot_max = alpha_f[ii][jj];}
  plot_min = floor(plot_min);
  plot_max = ceil(plot_max);
  /* min and max define the region of the axes in the PPM file 
     delta is the grid pitch: */
  delta = (plot_max-plot_min)/(pixel-1);
  sigma = delta/2.; /* Gauss parameter of the interpolation*/    /*P*/
  sigma2 = sigma * sigma;

  printf("plot_min = %4.2f deg; plot_max = %4.2f deg; \
delta = %4.4f deg; pixel = %i \n",
	 plot_min, plot_max, delta, pixel);

  for (i=1; i<=pixel; i++) {
    if(i%SCREEN_WRITE==0) printf("%i ",i); fflush(stdout);
    for (j=1; j<=pixel; j++) {
      sum = 0; norm = 0;
      for (ii=1; ii<=number; ii++) {
	di =  plot_min+(i-1)*delta - alpha_i[ii][1];
	if (fabs(di) < 3 * sigma) for (jj=0; jj<=1279; jj++) {    /*P*/
	  df =  plot_min+(j-1)*delta - alpha_f[ii][jj];
	  if (fabs(df) < 3 * sigma) {                             /*P*/
	    d2 = di * di + df * df;          
	    arg = -d2/sigma2;
	    if (arg > -LIM_EXP) {
	      w = exp(arg);
	      sum = sum + counts[ii][jj]*w;
	      norm = norm + w; } } } }
      PPM[i][j] = norm < EFFZERO ? 0 : sum/norm; 
      GNU[i][j] = PPM[i][j];
      a_i[i][j] = (plot_max*(i-1)+plot_min*(pixel-i))/(pixel-1);
      a_f[i][j] = (plot_max*(j-1)+plot_min*(pixel-j))/(pixel-1); } }
  printf("\n");
  return;}
/**********************    function: interpolate    ***************/


/**********************    function: logarithm    ***************/
void logarithm(){
  int i,j;
  
  for (i=1; i<=pixel; i++) 
    for (j=1; j<=pixel; j++)
      PPM[i][j]= log10(PPM[i][j] + BG);  
  return;}
/**********************    function: logarithm    ***************/


/****************************** ppm ********************************/
void ppm(void) { 
  FILE *file; 
  double dummy, fdummy, function_min, function_max;
  int i, j;

  /* against compiler warning message: */
   file = fopen(ppmfile, "w"); fclose( file);  

    /* calculation of the allover function range: */
    function_min = PPM[1][1]; 
    function_max = PPM[1][1]; 
    for (j=pixel; j>=1; j--)
      for (i=1; i<=pixel; i++) {
          if (function_min > PPM[i][j]) function_min = PPM[i][j];
          if (function_max < PPM[i][j]) function_max = PPM[i][j];}
    fprintf(stderr, "function range: min = %g; max = %g \n",
            function_min, function_max);
    if (fabs(function_max - function_min) < EFFZERO) 
      terminate("the function does not change within the range");
    
      file = fopen(ppmfile, "w");

      fprintf (file, "P3 \n");
      fprintf (file, "# ppm = portable pixmap \n");
      fprintf (file, "# pixel: R G B \n");
      fprintf (file, "# white : 255 255 255 \n");
      fprintf (file, "# black:    0   0   0 \n");
      fprintf (file, "# red:    255   0   0 \n");
      fprintf (file, "# green:    0 255   0 \n");
      fprintf (file, "# blue:     0   0 255  \n");
      fprintf (file, "#      black&white: P3->P1; pixel: D \n");
      fprintf (file, "# size: (e.g., 400 300) \n");
      fprintf (file, " %i   %i \n", pixel, pixel);
      fprintf (file, "# 8 bit/colour: \n");
      fprintf (file, "255 \n");
      fprintf (file, "# let us start: (e.g., 400x300=120 000 lines) \n");
      fprintf (file, "# (at the very end a newline character must be set)  \n");
      
      fdummy = (colour_max - colour_min)/(function_max-function_min);
      
      for (i = pixel; i >= 1; i--) 
        for (j = 1; j <= pixel; j++) {
	  dummy = (PPM[j][i] - function_min) * fdummy + colour_min; 
          colour (file, dummy, 1., 1.); }    
      /* closing with \n: */
      fprintf (file, " \n");
      fclose(file);

      message (ppmfile, "   has been written or rewritten.");

      sprintf(ppmfile, "%s%s", "scale", ".ppm");             /*P*/     
      file = fopen(ppmfile, "w");

      fprintf (file, "P3 \n");
      fprintf (file, "# ppm = portable pixmap \n");
      fprintf (file, "# pixel: R G B \n");
      fprintf (file, "# white : 255 255 255 \n");
      fprintf (file, "# black:    0   0   0 \n");
      fprintf (file, "# red:    255   0   0 \n");
      fprintf (file, "# green:    0 255   0 \n");
      fprintf (file, "# blue:     0   0 255  \n");
      fprintf (file, "#      black&white: P3->P1; pixel: D \n");
      fprintf (file, "# size: (e.g., 400 300) \n");
      fprintf (file, " %i   %i \n", pixel, pixel);
      fprintf (file, "# 8 bit/colour: \n");
      fprintf (file, "255 \n");
      fprintf (file, "# let us start: (e.g., 400x300=120 000 lines) \n");
      fprintf (file, "# (at the very end a newline character must be set)  \n");
      
      for (j = pixel; j >= 1; j--) 
	for (i = 1; i <= pixel; i++){
	  dummy = colour_min + 
	    (j-1.) * (colour_max - colour_min )/(pixel-1.); 
	  colour (file, dummy, 1., 1.); }
      /* closing with \n: */
      fprintf (file, " \n");
      fclose(file);
      

      message (ppmfile, "   has been written or rewritten.");
      message ("You may want to run the 'gqview' program",
	       "to display the 2D plot.\n");}
/*************************** end ppm *******************************/

/*************************    colour         ************ ********/
void colour (FILE *file, float h, float s, float v) {
  /* input -> RGB -> write into file 
     h: hue (0...1) colour family (this is the variable) (cyclic)
     s: saturation (0...1)
     v: value (0...1) (brightness)
     r,g,b: 0...1
     usual case:      h = variable, s=1; v=1;
     0.0 ... 0.25: red-yellow;  0.25 ... 0.5: green-blue
     0.5 ... 0.75 : blue;       0.75 ... 1.0: red */  
  int i;
  float f, p, q, t; 
  float r, g, b;  
  if (s==0) {/* achromatic - grey */ r=g=b=v; return;}
  h*=6; /* sector 0 to 5 */ i=floorf(h); 
  f=h-i; /* factorial part of h */
  p=v*(1-s); q=v*(1-s*f); t=v*(1-s*(1-f));
  switch (i){
  case 0: r=v; g=t; b=p; break;
  case 1: r=q; g=v; b=p; break;
  case 2: r=p; g=v; b=t; break;
  case 3: r=p; g=q; b=v; break;
  case 4: r=t; g=p; b=v; break;
  default:r=v; g=p; b=q; break;}

  fprintf (file, " %i  %i  %i  \n", (int)(floorf(255.*r)), 
           (int)(floorf(255.*g)), 
           (int)(floorf(255.*b))); }
/**************************** colour *******************************/


/****************************** specular ********************************/
double specular(double m, double alpha_i, double alpha_f) { 
  return fabs(sqrt(m*m+1)*(alpha_f-m*alpha_i));}
/************************    function: specular    *********************/


/****************************** horizontal ********************************/
double horizontal(double alpha_f0, double alpha_f) { 
  return fabs(alpha_f0-alpha_f);}
/************************    function: horizontal    *********************/


/****************************** vertical ********************************/
double vertical(double alpha_i0, double alpha_i) { 
  return fabs(alpha_i0-alpha_i);}
/************************    function: vertical    *********************/


/****************************** rocking ********************************/
double rocking(double alpha_i0, double alpha_f0, double alpha_i, double alpha_f) 
{return fabs((alpha_i0*alpha_f0-alpha_i0*alpha_f-alpha_f0*alpha_i)
	     /sqrt(alpha_i0*alpha_i0+alpha_f0*alpha_f0));}
/************************    function: rocking    *********************/


/***********************  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: alpha_f = m * alpha_i
                  distance: sqrt(1+m^2)*(alpha_f-m*alpha_i)
   all other:     equation: alpha_i/alpha_i0 + alpha_f/alpha_f0 = 1
                  distance:
 (alpha_i0^2+alpha_f0^2)*(alpha_i0*alpha_f0-alpha_i0*alpha_f-alpha_f0*alpha_i)
*/

void grafika(void){
  double data[MAX_PIXEL], width = 0.1, sum;                /*P*/
  char ch[MAX_LEN], dum, previous='s', present='s';
  double von, bis, von1, bis1, delta;  /* zoom parameters */
  int zoom=1;                          /* zoom factor */
  double zn=1;                         /* zoom range */ 
  FILE *file; 
  FILE *pipe;  
  char *kontr; 
  int i, j, logscale=1;
  double alpha_i0 = 0., alpha_f0= 0., m = 1, dummy;                        /*P*/
  
  sprintf(ch, "s");
  
  do {pipe = popen("gnuplot -geometry 900x400+10+10", "w");
    if (logscale == 1) fprintf (pipe, "set logscale y \n");
    else fprintf (pipe, "set nologscale y \n");
    fprintf (pipe, "set grid \n");
    fprintf (pipe, "set title '<%s>' \n", name);
    fflush(pipe);
    
  wrestart:
    
    switch (ch[0]){ 
	
    case '?': /* help */
      printf("\nRETURN lets the default value (displayed in [...]) unchanged \n"
	     "s: specular cut  slope = 1 \n" 
	     "r: general cut (intersections with the axes in degree) \n"
	     "h: horizontal cut (constant alpha_f in degree) \n"
	     "v: vertical   cut (constant alpha_i in degree) \n"
	     "l: logarithmic vertical axis in the 1D plot \n"
	     "n: linear      vertical axis in the 1D plot \n"
	     "p: makes postscript file \n"
	     "g: makes png        file \n"
	     "2: makes 2D gnuplot file \n"
	     "w: integration width in degree \n");
      break;

    case 'w': /* change the width */
      printf("integration width (degrees) [%8.4f] = ", width);
      fgets (ch, MAX_LEN, stdin);
      dummy=strtod(ch, &kontr);
      if (dummy!=0) width = dummy;
      ch[0] = previous;
      goto wrestart;
      break;

    case 'l': /* logscale on */
      logscale = 1;
      fprintf (pipe, "set logscale y \n");
      break;
      
    case 'n': /* logscale off */
      logscale = 0;
      fprintf (pipe, "set nologscale y \n");
      break;

    case '2': /* makes 2D gnuplot file */
      /* against compiler warning message: */
      file = fopen(gnuplot2Dfile, "w"); fclose( file);  
      file = fopen(gnuplot2Dfile, "w");
      for (i=1; i<=pixel; i++) {
	for (j=1; j<=pixel; j++) 
	  fprintf (file, " %f  \t %f \t %g \n", a_i[i][j], a_f[i][j], GNU[i][j]);
	  fprintf (file, " \n"); }
      fclose(file);
      printf ("\nThe file \"%s\" has been written or rewritten. \n",gnuplot2Dfile);
      break;

    case 'z':    /* zoom calculation: */
      sscanf (ch, "%1c%1i", &dum, &zoom); 
      if (zoom < 1) zoom = 1;
      zn = (zoom+1.)/2.;
      break;
      
    case '+':    /* zoom calculation: */
      zn = zn + .5; 
      break;
      
    case '=': /* zoom calculation: */
      zn = zn + .5; 
      break;
      
    case '-': /* zoom calculation: */
      zn = zn - .5;  
      if (zn < 1)   zn = 1;
      if (zn > zoom) zn = zoom;
      break;

    case 'p': /* output: postscript  */
      fprintf(pipe, "set terminal postscript color solid lw 2\n");      
      fprintf(pipe, "set output '%s' \n", psfile);  
      fflush(pipe);
      break;
      
    case 'g': /* output: png  */
      if (gnuplotversion == 0)
	{fprintf(pipe, "set terminal png \n");      
	  fprintf(pipe, "set output '%s' \n", pngfile);  
	  fflush(pipe);}
      break;
      
    case 's': /* specular */
      printf("slope [%8.4f] = ", m);
      fgets (ch, MAX_LEN, stdin);
      dummy = strtod(ch, &kontr);
      if (dummy!=0) m = dummy;
      for (i=1; i<=pixel; i++) {sum = 0;
	for (j=1; j<=pixel; j++)
	  if (specular (m,a_i[i][j],a_f[i][j]) < width) sum = sum + GNU[i][j];
	data[i] = sum; }
      file = fopen (gnuplotfile, "w");
      for (i=1; i<=pixel; i++)
	if (data[i]>PLOTZERO)	
	  fprintf(file, "%f  %f \n",a_i[i][1], data[i]);
      fclose(file);
      break;
      
    case 'v': /* vertical:   alpha_i=const. */                  
      printf("alpha_i0 [%8.4f] = ", alpha_i0);
      fgets (ch, MAX_LEN, stdin);
      dummy = strtod(ch, &kontr);
      if (dummy!=0) alpha_i0 = dummy;
      for (j=1; j<=pixel; j++) {sum = 0;
	for (i=1; i<=pixel; i++)
	  if (vertical (alpha_i0,a_i[i][j]) < width) sum = sum + GNU[i][j];
	data[j] = sum; }
      file = fopen (gnuplotfile, "w");
      for (i=1; i<=pixel; i++)
	if (data[i]>PLOTZERO)
	  fprintf(file, "%f  %f \n",a_f[1][i], data[i]);  
      fclose(file);
      break;
      
    case 'h': /* horizontal: alpha_f=const */  
      printf("alpha_f0 [%8.4f] = ", alpha_f0);
      fgets (ch, MAX_LEN, stdin);
      dummy = strtod(ch, &kontr);
      if (dummy!=0) alpha_f0 = dummy;
      for (i=1; i<=pixel; i++) {sum = 0;
	for (j=1; j<=pixel; j++)
	  if (horizontal (alpha_f0,a_f[i][j]) < width) sum = sum + GNU[i][j];
	data[i] = sum; }
      file = fopen (gnuplotfile, "w");
      for (i=1; i<=pixel; i++)
	if (data[i]>PLOTZERO)
	  fprintf(file, "%f  %f \n",a_i[i][1], data[i]);
      fclose(file);
      break;
      
    case 'r': /* rocking or similar */    
      printf("alpha_i0 [%8.4f] = ", alpha_i0);
      fgets (ch, MAX_LEN, stdin);
      dummy = strtod(ch, &kontr);
      if (dummy!=0) alpha_i0 = dummy;
      printf("alpha_f0 [%8.4f] = ", alpha_f0);
      fgets (ch, MAX_LEN, stdin);
      dummy = strtod(ch, &kontr);
      if (dummy!=0) alpha_f0 = dummy;
      for (i=1; i<=pixel; i++) {sum = 0;
	for (j=1; j<=pixel; j++)
	  if (rocking (alpha_i0,alpha_f0,a_i[i][j],a_f[i][j]) < width) 
	    sum = sum + GNU[i][j];
	data[i] = sum; }
      file = fopen (gnuplotfile, "w");
      for (i=1; i<=pixel; i++)
	if (data[i]>PLOTZERO)
	  fprintf(file, "%f  %f \n", a_i[i][pixel], data[i]);
      fclose(file);
      break;
      
    default: ;}
    
    von1 = alpha_i_min; 
    bis1 = alpha_i_max; 
    delta = (bis1 - von1)/zoom;
    von = von1 + (zn-1) * delta; bis = von + delta;
    
    /*    if (gnuplotversion == 0)
      fprintf (pipe,    
	       "plot [%f:%f][:] '%s' t '' w l lw 2 lt 1 \n",
	       von, bis, gnuplotfile);
    else
      fprintf (pipe,   
	       "plot [%f:%f][:] '%s' t '' w l lw 2    1 \n",
	       von, bis, gnuplotfile);
	       fflush(pipe); */
    
    if (gnuplotversion == 0)
      fprintf (pipe,    
	       "plot [:][:] '%s' t '' w l lw 2 lt 1 \n", gnuplotfile);
    else
      fprintf (pipe,   
	       "plot [:][:] '%s' t '' w l lw 2    1 \n", gnuplotfile);
    fflush(pipe); 

    if (ch[0] == 'p') 
      printf ("\nThe file \"%s\" has been written or rewritten. \n", psfile);
      
    if (ch[0] == 'g') 
      printf ("\nThe file \"%s\" has been written or rewritten. \n", pngfile);
          
    printf
      ("\nlogplot: 'l'  linplot: 'n'  integration width in degrees [%5.3f]: 'w'\n" 
       "makes ps file: 'p'  makes png file: 'g'  makes 2D gnuplot file: '2'\n"
       "help: '?'  quit: 'q'\n"
       "specular: 's'  horizontal: 'h'  vertical: 'v'  rocking: 'r' --->", width);
    

    fgets (ch, MAX_LEN, stdin);      
    if (present!='w') previous = present;
    present = ch[0];   

    fprintf (pipe, "q  \n");        fflush(pipe);
  } while (ch[0] != 'q');
  fclose (pipe);}
/*************************  grafika ***************************/


/**********************    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    ***************/


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


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