#define DATUM                  "2009_06_29"
/****************************************************
This is the file 'mythen.c.2009_06_29', 
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) at the 
"Large Four Circle 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:    
   cc -lm mythen.c -Wall -o ~/bin/mythen

   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_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);

/* global declarations: */
int    counts [MAX_FILE][MAX_CHANNEL];  /* experimental results */
double q_x[MAX_FILE][MAX_CHANNEL];
double q_y[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 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 *spec_file;            /* 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];
int scan_Nr;

int m_all, m_act; /*  number of files: mythen_all, mythen_actual*/

/* 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;
  message ("\nProgram version:   mythen.c." DATUM, "" );
  if (argc != 3) {message("\n You have made a mistake in the input,",
                          "please try again using the format:\n");
    message (" ", 
 "mythen <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  ");}
  
  spec_file = argv[1];
  message("filename=", spec_file);
  
  scan_Nr = strtod(argv[2], &kontr);
  if(kontr == argv[2]) terminate(" wrong second parameter\n");
  
  if (scan_Nr >= MAX_FILE){
    printf("The number of files >= MAX_FILE= %i:\n", MAX_FILE);
    terminate("please correct the C code.");}

   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 */  
  read_spec_file();
   read_mythen_files ();
   q_calculation ();  
   write_gnuplot2D();
   grafika();
   message ("\n\nProgram version:   mythen.c." DATUM, "\n" );
 
   error ("terminal"); 
  return 0;}
/**********************    function: main    ***************/



/**********************    function: read_spec_file    ***************/
void read_spec_file () {
  /* */
  //  char filename [MAX_LEN];
  char line [MAX_LEN+1];
  char dummy [MAX_LEN];
  char motor [MAX_LEN] ;
  float start,stop,step, time;  
 char motor1 [MAX_LEN],motor2[MAX_LEN];
 float start1,start2,stop1,stop2;  
  int idummy ; 
  FILE *file; 

  //  sprintf(filename, name);
  // printf("read_file_anfang: filename=%s \n", filename);

  
file=fopen(spec_file,"r");
  if (file==NULL)  terminate(" The desired file does not exist.  ");
  while (fgets (line, MAX_LEN, file) != NULL) {
    if ((line[0]=='#')&&(line[1]=='S')){
      sscanf (line, "%s%i%s", dummy,&idummy,dummy);
      //printf("scan=%i    scanA=%i  dummy1=%s dummy2=%s \n",scan,scanA,&dummy1,&dummy2);
      if (idummy==scan_Nr){
	printf("scanNr.=: %i     scan=:  %s\n",scan_Nr,dummy);
	if (strcmp(dummy, "ascan")); else {
	  sscanf (line, "%s%s%s%s%f%f%f%f", dummy, dummy, dummy,motor,&start,&stop,&step,&time);
	  m_all=step+1;
	printf("m_all=%i:\n",m_all) ;
	}
	if (strcmp(dummy, "a2scan")); else {
	  sscanf (line, "%s%s%s%s%f%f%s%f%f%f%f", dummy, dummy, dummy,motor1,&start1,&stop1,motor2,&start2,&stop2,&step,&time);
	  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%s", dummy, dummy, dummy,dummy, dummy, dummy,dummy, dummy, dummy,&step,dummy);
	  m_all=step+1;
	printf("m_all=%i:\n",m_all) ;
	}

	break;     }   } }

  while (fgets (line, MAX_LEN, file) != NULL) {
    if ((line[0]=='#')&&(line[1]=='L')){
     
      //sscanf (line, "%s%s%s", dummy,dummy,dummy);
      //printf("L-line=: %s   \n",dummy);
      //printf("scan=%i    scanA=%i  dummy1=%s dummy2=%s \n",scan,scanA,&dummy1,&dummy2);
      //
	break;     }   } 
 m_act=0;
 while (fgets (line, MAX_LEN, file) != NULL) {
   if ((line[0]=='#')||(strlen(line)==1))  break;
  m_act=m_act+1;  
  } 

 printf("m_act=: %i\n",m_act);


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


/**********************    function: read_mythen_files    ***************/
void read_mythen_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<=m_act; ii++){
    sprintf(filename, "%s%s%s%i%s%i%s%i%s", spec_file,"_", "scan",scan_Nr,"_",ii,"of", m_all, ".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_mythen_files    ***************/

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

  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<=1279; jj++){                           
	q_x[ii][jj] = ii;
	q_y[ii][jj] = jj;}          /*   or	q_y[ii][jj] = 1279-jj     */   
    break;

   case 1:                                    /* tth = th */
    /* sample horizon: the same channel */
    for (ii=1; ii<=m_act; ii++)
      for (jj=0; jj<=1279; 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 */
    /* specular direction: the same channel */
    for (ii=1; ii<=m_act; ii++)
      for (jj=0; jj<=1279; 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;
    
  default: terminate ("geometry parameter (GEO): illegal value ");
    break;}
  
  return;}
/**********************    function: q_calculation    ***************/



/* **********************  write_gnuplot2D******************/

void  write_gnuplot2D(void){
  FILE *file; 
  int i,j;
 /* makes 2D gnuplot file */
    /* against compiler warning message: */
  printf("m_act=%i \n  ",m_act);
  file = fopen(gnuplot2Dfile, "w"); fclose( file);  
  file = fopen (gnuplot2Dfile, "w") ;
  for (i=1; i<=m_act; i++){ 
    for (j=0; j<=1279; j++) 
      fprintf (file, " %f \t  %f \t  %i  \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",gnuplot2Dfile);
}
/* **********************  write_gnuplot2D******************/


/***********************  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];
  FILE *pipe;  
  int logscale=1;  
  /* gnuplot settings */  
  pipe = popen("gnuplot -geometry 900x900+10+10", "w");
 start:
  if (logscale == 1) {fprintf (pipe, "set logscale z \n");
    fprintf (pipe, "set logscale cb \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");
  fprintf (pipe, "set title '<%s>' \n", spec_file);
  fflush(pipe);
 

    /* gnuplot command */
  fprintf (pipe,    
	   "splot [:][:][:] '%s' t '' with pm3d \n", gnuplot2Dfile);
  fflush(pipe); 
  fprintf(stderr, "'li' = linear scale in z;  'lo' = logarithmic scale; 'q' = QUIT; Further with RETURN: "); 
  fgets(ch,MAX_LEN,stdin);
  if ((ch[0]=='l') &&(ch[1]=='i')) {logscale=0; goto start;} else
    if ((ch[0]=='l') &&(ch[1]=='o')) {logscale=1; goto start;} else
    if ((ch[0]!='q')&&(ch[0]!='Q')) goto start; 
  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     ********************/