#define DATUM "2010_06_07" /**************************************************** This is the file 'file_read.c.2010_06_07', genuine name: 'file_read.c'. The aim of this software is the analysis of the experimental results obtained at the beamtimes in July 2009 December 2009 at the surface diffraction beamline of the Max-Planck-Institut fuer Metallforschung at the ANKA synchrotron facility in Karlsruhe. The program 'file_read' was written at the Max-Planck-Institut fuer Metallforschung (Stuttgart). Authors: Agnes Szoekefalvi-Nagy, Janos Major (major@mf.mpg.de) The newest version of the source code 'file_read.c' can freely be obtained from the web page of the Max Planck Institute for Metals Research (www.mf.mpg.de), page of the Department for Low Dimensional and Metastable Materials, page "software". The program 'file_read' 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. The included Fortran fit routine minuit is a GPL (GNU General Public Licence, http://www.gnu.org/copyleft/gpl.html) software of CERN (Geneva, Switzerland) (http://consult.cern.ch/writeup/minuit/). The c-minuit coupling of the Fortran minuit routine with the C language program, the c-minuit system was developed by Andras Major (http://c-minuit.sourceforge.net). For installation of the minuit C bibliothek see the man page of the superfit program and the file c-minuit-CVS-20031201.tar.gz (web page of the Max Planck Institute for Metals Research (www.mf.mpg.de), page of the Department for Low Dimensional and Metastable Materials, page "software"). 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: file_read.c.2009_08_14 The standard site for this file is: $HOME/c-minuit// To compile and link please use: make install Adjustable parameters in the lines with: */ /*P*/ /* *************************************************/ /* for monitoring all errors: int ERROR = 0; for stopping at the first error: int ERROR = 1; */ int ERROR = 0; #define __NO_MATH_INLINES /* makes slower but, may be, more precise */ #include #include #include #include #include #include #include "cfortran.h" #include "minuitfcn.h" #include "minuit.h" #define EICHUNG 0.00804 /* mythen detector: degree/channel */ /*P*/ #define MAX_FILE 5100 /* max. number of files */ #define MAX_PIXEL 5000 /* 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 NEARZERO 1e-10 /* zero for the matrix functions */ #define LIM_EXP 70. /* limit for the function exp */ #define SCREEN_WRITE 50 /* 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 #define MAX_MESSPUNKT 10000000 /* max. number of measuring points */ #define MAX_FIT 5000 /* max. number of fitted plots */ #define NM 1000 /* max. number of fit parameters */ #define MODEL model3 /* which model function will bw used */ /*P*/ /* function prototypes: */ void message (char *text, char *text1); void error (char *text); void terminate (char *text); void wait (char *text); void read_dat_files(void); void read_fit_dat_files(void); void gnuplot_read_dat_files_t (void); void gnuplot_read_dat_files_q (void); double model1(double); double model2(double); double model3(double); void fcn(int npar,double* grad,double* fcnval,double* x,int iflag,void* futil); int tr(int); double power(double, double); void header (void); /* global declarations: */ int mode; double start[NM], step[NM], min[NM], max[NM]; /* do not change this line */ char parname[NM][MAX_LEN]; /* names of the MINUIT pars */ double fcn_end, fcnred_end; /* minfit results */ int n_par; /* minfit results */ double HH,KK,LL,MON[5000],DET[5000], MMON[5000], MMONITOR; /* fit data */ double S1, S2, S3, S4; double A1_1, A1_2, A1_3, A1_4, T1_1, T1_2, T1_3, T1_4; double A2_1, A2_2, A2_3, A2_4, T2_1, T2_2, T2_3, T2_4; double A3_1, A3_2, A3_3, A3_4, T3_1, T3_2, T3_3, T3_4; double D3_1, D3_2, D3_3, D3_4; double S1f[MAX_FIT], S2f[MAX_FIT], S3f[MAX_FIT], S4f[MAX_FIT]; double A1_1f[MAX_FIT], A1_2f[MAX_FIT], A1_3f[MAX_FIT], A1_4f[MAX_FIT]; double T1_1f[MAX_FIT], T1_2f[MAX_FIT], T1_3f[MAX_FIT], T1_4f[MAX_FIT]; double A2_1f[MAX_FIT], A2_2f[MAX_FIT], A2_3f[MAX_FIT], A2_4f[MAX_FIT]; double T2_1f[MAX_FIT], T2_2f[MAX_FIT], T2_3f[MAX_FIT], T2_4f[MAX_FIT]; double A3_1f[MAX_FIT], A3_2f[MAX_FIT], A3_3f[MAX_FIT], A3_4f[MAX_FIT]; double T3_1f[MAX_FIT], T3_2f[MAX_FIT], T3_3f[MAX_FIT], T3_4f[MAX_FIT]; double D3_1f[MAX_FIT], D3_2f[MAX_FIT], D3_3f[MAX_FIT], D3_4f[MAX_FIT]; double Hf [MAX_FIT], Kf [MAX_FIT], Lf [MAX_FIT]; double t0f[MAX_FIT], t1f[MAX_FIT], t2f[MAX_FIT], t3f[MAX_FIT], t4f[MAX_FIT]; double FCNRED[MAX_FIT]; 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 */ double H[MAX_MESSPUNKT], K[MAX_MESSPUNKT], L[MAX_MESSPUNKT]; double maxdm=0; /* MAX(det/mon) */ int Time[MAX_MESSPUNKT], mon[MAX_MESSPUNKT], det[MAX_MESSPUNKT]; char *datfile; /* first (common) part of the file names */ char *datpath; /* path 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 first, last, total; double t0,t1,t2,t3,t4; /* times of voltage chages: fitted values */ double t0n, t1n, t2n,t3n,t4n; /* times of voltage chages: nominal values */ double COUNTSperMMON [5000][5000]; /* exp. counts[shot][time]*/ int errorcount = 0; /* for the function error */ int gnuplotversion; /* 0 if >= 4.0; >0 if <=3.7 */ int fcncount=0; char headerchar[2000]; /********************** function: main ***************/ int main(int argc, char *argv[]){ char *kontr; if (argc != 12) {message("\n You have made a mistake in the input,", "please try again using the format:\n"); message (" ", "file_read \n"); message (" =","1: reads data files and makes nothing else"); message (" ","2: makes 1D graphic (png) files by gnuplot "); message (" "," (axes: t, counts/mean_monitor) "); message (" "," (for autoscaling please edit the source code)"); message (" ","3: makes 2D graphic (png) files by gnuplot "); message (" "," (axes: H, K; colour: counts/mean_monitor) "); message (" "," (for autoscaling please edit the source code)"); message (" ","4: reads data files and fits by minuit "); message (" "," (counts/mean_monitor values are fitted; "); message (" "," mean_monitor = average within the shot) "); message (" "," 1-: all files considered "); message (" "," gnuplot2D files will be made "); message (" "," (axes: H, K; colour: fit parameters) "); message (" "," else: on screen graphical test of the fits "); message (" =","path of the data files (e.g., . or ../ or DAT)"); message (" =","first (common) part of the file names"); message (" ="," 0 ... t0n: voltage = 0 (t0n: nominal value)"); message (" =","t0n ... t1n: voltage > 0 (t1n: nominal value)"); message (" =","t1n ... t2n: voltage = 0 (t2n: nominal value)"); message (" =","t2n ... t3n: voltage < 0 (t3n: nominal value)"); message (" =","t3n ... t4n: voltage = 0 (t4n: nominal value)"); message (" ","t4n is the periode time of the voltage sequence"); message (" =","number of mesh points"); message (" =","first data file to be considered"); message (" =","last data file to be considered"); message ( "Attention: all files produced by the present program will be located in\n", " the directory where the program is started. Existing files with"); message ( " the same name will be overwritten without warning.",""); fprintf (stderr, "Program version: file.read." DATUM "\n"); terminate (" input error ");} mode = strtod(argv[1], &kontr); if(kontr == argv[1]) terminate(" wrong first parameter\n"); datpath = argv[2]; message("path of the data files= ", datpath); datfile = argv[3]; message("common part of the file names=", datfile); t0n = strtod(argv[4], &kontr); if(kontr == argv[4]) terminate(" wrong 4th parameter\n"); t1n = strtod(argv[5], &kontr); if(kontr == argv[5]) terminate(" wrong 5th parameter\n"); t2n = strtod(argv[6], &kontr); if(kontr == argv[6]) terminate(" wrong 6th parameter\n"); t3n = strtod(argv[7], &kontr); if(kontr == argv[7]) terminate(" wrong 7th parameter\n"); t4n = strtod(argv[8], &kontr); if(kontr == argv[8]) terminate(" wrong 8th parameter\n"); total = strtod(argv[9], &kontr); if(kontr == argv[9]) terminate(" wrong 9th parameter\n"); first = strtod(argv[10], &kontr); if(kontr == argv[10]) terminate(" wrong 10th parameter\n"); last = strtod(argv[11], &kontr); if(kontr == argv[11]) terminate(" wrong 11th parameter\n"); t0 = t0n; t1 = t1n; t2 = t2n; t3 = t3n; t4 = t4n; /* 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 (mode){ case 1: /* this case if only for a general test */ read_dat_files (); break; case 2: read_dat_files (); /* maxdm (range of the counts/mon) and MMON[ii] = mean_monitor will be defined */ gnuplot_read_dat_files_t (); break; case 3: read_dat_files (); /* maxdm (range of the counts) and MMON[ii] = mean_monitor will be defined */ gnuplot_read_dat_files_q (); break; case 4: read_fit_dat_files (); break; default: terminate(" Illegal mode value. "); break;} message ("Program version: file_read.c." DATUM, "\n" ); error ("terminal"); return 0;} /********************** function: main ***************/ /*************************** header ****************************/ void header (void) {/* prepares the string headerfile to be copied as a header for the out files */ char c[2000]; time_t curtime; curtime = time (NULL); strcpy(c, "# "); strcat(headerchar,c); strcpy(c, asctime (localtime (&curtime))); strcat(headerchar,c); strcpy(c, "# First (common) part of the data file names: \n"); strcat(headerchar,c); strcpy(c, "# "); strcat(headerchar,c); strcpy(c, datpath); strcat(headerchar,c); strcpy(c, datfile); strcat(headerchar,c); strcpy(c, "\n"); strcat(headerchar,c); sprintf (c, "# The initial time values in experimental units: \n"); strcat(headerchar,c); sprintf (c, "# t0n=%6.2f t1n=%6.2f t2n=%6.2f t3n=%6.2f t4n=%6.2f \n", t0n, t1n, t2n, t3n, t4n); strcat(headerchar,c); sprintf (c, "# Further initial parameters: \n"); strcat(headerchar,c); sprintf (c, "# mode=%1i total=%4i first=%4i last=%4i \n", mode, total, first, last); strcat(headerchar,c); strcpy(c, "# This file was produced by the software 'file-read' version " DATUM "\n"); strcat(headerchar,c); if (strcmp("MODEL","model1")); strcpy(c, "# The used model function for the fit is the 'model1',\n" "# which is the stretched exponential fuction \n"); if (strcmp("MODEL","model2")); strcpy(c, "The used model function for the fit is the 'model2',\n" "# which is the sum of two exponentials \n"); if (strcmp("MODEL","model3")); strcpy(c, "# The used model function for the fit is the 'model3',\n" "# which is the sum of three exponentials;\n" "# the third one can be delayed and gives zero for negative \n" "# arguments (Heaviside function) \n"); strcat(headerchar,c); strcpy(c, "# \n"); strcat(headerchar,c);} /*************************** header ****************************/ /*************************** model1 ****************************/ double model1(double t) { /* model function: stretched exponential */ double dummy; if (t<=t0) {dummy = -power((t+t4-t3)/T1_4, A2_4); if (dummy < -100) dummy =S4; else dummy = S4*(1. + A1_4*exp(dummy));} else if (t<=t1) {dummy = -power( (t-t0)/T1_1, A2_1); if (dummy < -100) dummy =S1; else dummy = S1*(1. + A1_1*exp(dummy));} else if (t<=t2) {dummy = -power( (t-t1)/T1_2, A2_2); if (dummy < -100) dummy =S2; else dummy = S2*(1. + A1_2*exp(dummy));} else if (t<=t3) {dummy = -power( (t-t2)/T1_3, A2_3); if (dummy < -100) dummy =S3; else dummy = S3*(1. + A1_3*exp(dummy));} else {dummy = -power( (t-t3)/T1_4, A2_4); if (dummy < -100) dummy =S4; else dummy = S4*(1. + A1_4*exp(dummy));} /* printf ("t= %g model = %g \n",t, dummy); */ return (dummy);} /*************************** model1 ****************************/ /*************************** model2 ****************************/ double model2(double t) { /* model function: sum of two exponentials */ double dummy, Dummy; if (t -100) dummy = dummy + S4*A2_4*exp(Dummy); } else if (t -100) dummy = dummy + S1*A2_1*exp(Dummy); } else if (t -100) dummy = dummy + S2*A2_2*exp(Dummy); } else if (t -100) dummy = dummy + S3*A2_3*exp(Dummy); } else { Dummy = -(t-t3)/T1_4; if (Dummy < -100) dummy =S4; else dummy = S4*(1. + A1_4*exp(Dummy)); Dummy = -(t-t3)/T2_4; if (Dummy > -100) dummy = dummy + S4*A2_4*exp(Dummy); } /* printf ("t= %g model = %g \n",t, dummy); */ return (dummy);} /*************************** model2 ****************************/ /*************************** model3 ****************************/ double model3(double t) { /* model function: - sum of three exponentials - the third one can be delayed and gives zero for negative arguments (Heaviside function) */ double dummy, Dummy; if (t -100) dummy = dummy + S4*A1_4*exp(Dummy); Dummy = -(t+t4-t3)/T2_4; if (Dummy > -100) dummy = dummy + S4*A2_4*exp(Dummy); Dummy = -(t+t4-t3-D3_4)/T3_4; if ((Dummy > -100)&&(Dummy < 0)) dummy = dummy + S4*A3_4*exp(Dummy); } else if (t -100) dummy = dummy + S1*A1_1*exp(Dummy); Dummy = -(t-t0)/T2_1; if (Dummy > -100) dummy = dummy + S1*A2_1*exp(Dummy); Dummy = -(t-t0-D3_1)/T3_1; if ((Dummy > -100)&&(Dummy < 0)) dummy = dummy + S1*A3_1*exp(Dummy); } else if (t -100) dummy = dummy + S2*A1_2*exp(Dummy); Dummy = -(t-t1)/T2_2; if (Dummy > -100) dummy = dummy + S2*A2_2*exp(Dummy); Dummy = -(t-t1-D3_2)/T3_2; if ((Dummy > -100)&&(Dummy < 0)) dummy = dummy + S2*A3_2*exp(Dummy); } else if (t -100) dummy = dummy + S3*A1_3*exp(Dummy); Dummy = -(t-t2)/T2_3; if (Dummy > -100) dummy = dummy + S3*A2_3*exp(Dummy); Dummy = -(t-t2-D3_3)/T3_3; if ((Dummy > -100)&&(Dummy < 0)) dummy = dummy + S3*A3_3*exp(Dummy); } else { dummy =S4; Dummy = -(t-t3)/T1_4; if (Dummy > -100) dummy = dummy + S4*A1_4*exp(Dummy); Dummy = -(t-t3)/T2_4; if (Dummy > -100) dummy = dummy + S4*A2_4*exp(Dummy); Dummy = -(t-t3-D3_4)/T3_4; if ((Dummy > -100)&&(Dummy < 0)) dummy = dummy + S4*A3_4*exp(Dummy); } /* printf ("t= %g model = %g \n",t, dummy); */ return (dummy);} /*************************** model3 ****************************/ /**************************** fcn *****************************/ void fcn(int npar,double* grad,double* fcnval,double* x,int iflag,void* futil) { int i; S1 = x[ 0]; S2 = x[ 1]; S3 = x[ 2]; S4 = x[ 3]; A1_1 = x[ 4]; A1_2 = x[ 5]; A1_3 = x[ 6]; A1_4 = x[ 7]; T1_1 = x[ 8]; T1_2 = x[ 9]; T1_3 = x[10]; T1_4 = x[11]; A2_1 = x[12]; A2_2 = x[13]; A2_3 = x[14]; A2_4 = x[15]; T2_1 = x[16]; T2_2 = x[17]; T2_3 = x[18]; T2_4 = x[19]; A3_1 = x[20]; A3_2 = x[21]; A3_3 = x[22]; A3_4 = x[23]; T3_1 = x[24]; T3_2 = x[25]; T3_3 = x[26]; T3_4 = x[27]; D3_1 = x[28]; D3_2 = x[29]; D3_3 = x[30]; D3_4 = x[31]; t0 = x[32]; t1 = x[33]; t2 = x[34]; t3 = x[35]; t4 = x[36]; /* here begins the function to be minimized (*fcnval=) : */ *fcnval=0.; /* Maximum Likelihood fit with Poisson statistics */ for (i=1; i<=floor(t4); i++) {*fcnval+=MODEL(i)*MMONITOR-DET[i]; if (DET[i]>1.e-6) *fcnval+=DET[i]*log(DET[i]/MODEL(i)/MMONITOR); } *fcnval = *fcnval*2.; /* print to screen: fit progress: fcncount = fcncount + 1; if (fcncount%SCREEN_WRITE==0) printf("fcn status: (%4i) fcn = %12.4f fcn_red = %12.4f \n", fcncount, *fcnval, *fcnval/(t4-npar)); */ if (iflag==3){fcn_end = *fcnval; fcnred_end = *fcnval/(t4-npar); n_par=npar;}} /**************************** fcn *****************************/ /********************** function: read_fit_dat_files ***************/ void read_fit_dat_files () { /* experimental counts will be read in: */ int ii, mtime, lineindex, i, j; char filename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; /* segedvaltozo */ int imon, idet; /* segedvaltozo egesztipusu */ char Dummy[MAX_LEN]; int dummy; double ddummy; int idummy, maxcall, parameter; char strategy[20], sdummy[MAX_LEN+1]; char ch[MAX_LEN]; int fitnumber; FILE *file; FILE *pipe; for (ii=first; ii<=last; ii++) { fcncount = 0; sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); /* printf("%s \n", filename); */ printf("%i ", ii); fflush(stdout); file=fopen(filename,"r"); mtime=0; lineindex=0; MMONITOR = 0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) sscanf (line, "%s%f%f%f", Dummy, &fH, &fK, &fL); else if (lineindex>12) {sscanf (line, "%i%i", &imon, &idet); mtime=mtime+1; HH=fH; KK=fK; LL=fL; MON[mtime]=imon; MMONITOR = MMONITOR + imon; DET[mtime]=idet; } } MMONITOR = MMONITOR/t4; MMON[ii] = MMONITOR; fclose(file); MNINIT(5,6,7); /* do not change this line, please */ /* t0 = t0n; t1 = t1n; t2 = t2n; t3 = t3n; t4 = t4n; = times in microseconds: nominal switching times */ sprintf(parname[ 1],"S1");start[ 1]=1;step[ 1]=.1;min[ 1]=0;max[ 1]=100; sprintf(parname[ 2],"S2");start[ 2]=1;step[ 2]=.1;min[ 2]=0;max[ 2]=100; sprintf(parname[ 3],"S3");start[ 3]=1;step[ 3]=.1;min[ 3]=0;max[ 3]=100; sprintf(parname[ 4],"S4");start[ 4]=1;step[ 4]=.1;min[ 4]=0;max[ 4]=100; sprintf(parname[ 5],"A1_1");start[ 5]=0.;step[ 5]=.01;min[ 5]=-.5;max[ 5]=5; sprintf(parname[ 6],"A1_2");start[ 6]=0.;step[ 6]=.01;min[ 6]=-.5;max[ 6]=5; sprintf(parname[ 7],"A1_3");start[ 7]=0.;step[ 7]=.01;min[ 7]=-.5;max[ 7]=5; sprintf(parname[ 8],"A1_4");start[ 8]=0.;step[ 8]=.01;min[ 8]=-.5;max[ 8]=5; sprintf(parname[ 9],"T1_1");start[ 9]=5;step[ 9]=.1;min[ 9]=.01;max[ 9]=5; sprintf(parname[10],"T1_2");start[10]=5;step[10]=.1;min[10]=.01;max[10]=5; sprintf(parname[11],"T1_3");start[11]=5;step[11]=.1;min[11]=.01;max[11]=5; sprintf(parname[12],"T1_4");start[12]=5;step[12]=.1;min[12]=.01;max[12]=5; sprintf(parname[13],"A2_1");start[13]=.0;step[13]=.1;min[13]=-.5;max[13]=5; sprintf(parname[14],"A2_2");start[14]=.0;step[14]=.1;min[14]=-.5;max[14]=5; sprintf(parname[15],"A2_3");start[15]=.0;step[15]=.1;min[15]=-.5;max[15]=5; sprintf(parname[16],"A2_4");start[16]=.0;step[16]=.1;min[16]=-.5;max[16]=5; sprintf(parname[17],"T2_1");start[17]=75;step[17]=.1;min[17]=10.;max[17]=4000; sprintf(parname[18],"T2_2");start[18]=75;step[18]=.1;min[18]=10.;max[18]=4000; sprintf(parname[19],"T2_3");start[19]=75;step[19]=.1;min[19]=10.;max[19]=4000; sprintf(parname[20],"T2_4");start[20]=75;step[20]=.1;min[20]=10.;max[20]=4000; sprintf(parname[21],"A3_1");start[21]=.0;step[21]=.0;min[21]=-.0;max[21]=.5; sprintf(parname[22],"A3_2");start[22]=.0;step[22]=.0;min[22]=-.0;max[22]=.5; sprintf(parname[23],"A3_3");start[23]=.0;step[23]=.0;min[23]=-.0;max[23]=.5; sprintf(parname[24],"A3_4");start[24]=.0;step[24]=.0;min[24]=-.0;max[24]=.5; sprintf(parname[25],"T3_1");start[25]=50;step[25]=.0;min[25]=.01;max[25]=2; sprintf(parname[26],"T3_2");start[26]=50;step[26]=.0;min[26]=.01;max[26]=2; sprintf(parname[27],"T3_3");start[27]=50;step[27]=.0;min[27]=.01;max[27]=2; sprintf(parname[28],"T3_4");start[28]=50;step[28]=.0;min[28]=.01;max[28]=2; sprintf(parname[29],"D3_1");start[29]=0;step[29]=.0;min[29]=1;max[29]=200; sprintf(parname[30],"D3_2");start[30]=0;step[30]=.0;min[30]=1;max[30]=200; sprintf(parname[31],"D3_3");start[31]=0;step[31]=.0;min[31]=1;max[31]=200; sprintf(parname[32],"D3_4");start[32]=0;step[32]=.0;min[32]=1;max[32]=200; sprintf(parname[33],"t0");start[33]=t0n;step[33]=0;min[33]=t0n-5;max[33]=t0n+5; sprintf(parname[34],"t1");start[34]=t1n;step[34]=0;min[34]=t1n-5;max[34]=t1n+5; sprintf(parname[35],"t2");start[35]=t2n;step[35]=0;min[35]=t2n-5;max[35]=t2n+5; sprintf(parname[36],"t3");start[36]=t3n;step[36]=0;min[36]=t3n-5;max[36]=t3n+5; sprintf(parname[37],"t4");start[37]=t4n;step[37]=0;min[37]=t4n-5;max[37]=t4n+5; for (j=1; j<=37; j++) MNPARM(j, parname[j], start[j], step[j], min[j], max[j], dummy); parameter = 2; /* 0: minimum number of fcn calls 1: (minuit default) 2: most precise results */ sprintf(strategy, "%s %i", "SET STRATEGY", parameter); MNCOMD(minuitfcn, strategy, dummy, NULL); /* default: SET WARNINGS SET NOWARNINGS */ sprintf(strategy, "%s", "SET NOWARNINGS"); MNCOMD(minuitfcn, strategy, dummy, NULL); parameter = 0; /* 0: minimum output 1: default value, normal output 2: additional output 3: maximum output*/ sprintf(strategy, "%s %i", "SET PRINTOUT", parameter); MNCOMD(minuitfcn, strategy, dummy, NULL); parameter = 1; /* 0.5: -log Likelihood (1 sigma) 1 : (minuit default) chisquare function (1 sigma) 4 : chisquare function (2 sigma) */ sprintf(strategy, "%s %i", "SET ERRORDEF", parameter); MNCOMD(minuitfcn, strategy, dummy, NULL); maxcall = 5000; sprintf(strategy, "%s %i", "SIMPLEX", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ maxcall = 5000; sprintf(strategy, "%s %i", "MINIMIZE", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ sprintf(strategy, "%s %i", "CALL", 3); MNCOMD(minuitfcn, strategy, dummy, NULL); /* results of the fit */ for (j=1; j<=25; j++) /* fit results */ MNPOUT(j, sdummy, start[j], step[j], ddummy, ddummy, idummy); S1f[ii] = start[ 1]; S2f[ii] = start[ 2]; S3f[ii] = start[ 3]; S4f[ii] = start[ 4]; A1_1f[ii] = start[ 5]; A1_2f[ii] = start[ 6]; A1_3f[ii] = start[ 7]; A1_4f[ii] = start[ 8]; T1_1f[ii] = start[ 9]; T1_2f[ii] = start[10]; T1_3f[ii] = start[11]; T1_4f[ii] = start[12]; A2_1f[ii] = start[13]; A2_2f[ii] = start[14]; A2_3f[ii] = start[15]; A2_4f[ii] = start[16]; T2_1f[ii] = start[17]; T2_2f[ii] = start[18]; T2_3f[ii] = start[19]; T2_4f[ii] = start[20]; A3_1f[ii] = start[21]; A3_2f[ii] = start[22]; A3_3f[ii] = start[23]; A3_4f[ii] = start[24]; T3_1f[ii] = start[25]; T3_2f[ii] = start[26]; T3_3f[ii] = start[27]; T3_4f[ii] = start[28]; D3_1f[ii] = start[29]; D3_2f[ii] = start[30]; D3_3f[ii] = start[31]; D3_4f[ii] = start[32]; t0f [ii] = start[33]; t1f [ii] = start[34]; t2f [ii] = start[35]; t3f [ii] = start[36]; t4f [ii] = start[37]; Hf[ii] = HH; Kf[ii] = KK; Lf[ii] = LL; FCNRED[ii] = fcnred_end; printf("ii=%i n_par=%i fcnred=%g \n",ii,n_par,fcnred_end);} if ((first > 1)||(last < total)) { /* Compare the data and the fit results by gnuplot only if not the complete data set is compared: */ fitnumber = first; pipe = popen("gnuplot -geometry 900x900+10+10", "w"); fprintf (pipe, "set term wxt noraise \n"); do { printf("number for comparison [range: %i-%i; quit: 'x']: ", first, last); fgets (ch, MAX_LEN, stdin); if (ch[0]=='x') return; if ((ch[0]=='+')||(ch[0]=='=')||(ch[0]==' ')) fitnumber=fitnumber+1; else if ((ch[0]=='-')||(ch[0]=='_')) fitnumber=fitnumber-1; else sscanf (ch, "%i", &fitnumber); if (fitnumber < first) fitnumber = first; if (fitnumber > last) fitnumber = last; S1 = S1f [fitnumber]; S2 = S2f [fitnumber]; S3 = S3f [fitnumber]; S4 = S4f [fitnumber]; A1_1 = A1_1f [fitnumber]; A1_2 = A1_2f [fitnumber]; A1_3 = A1_3f [fitnumber]; A1_4 = A1_4f [fitnumber]; T1_1 = T1_1f [fitnumber]; T1_2 = T1_2f [fitnumber]; T1_3 = T1_3f [fitnumber]; T1_4 = T1_4f [fitnumber]; A2_1 = A2_1f [fitnumber]; A2_2 = A2_2f [fitnumber]; A2_3 = A2_3f [fitnumber]; A2_4 = A2_4f [fitnumber]; T2_1 = T2_1f [fitnumber]; T2_2 = T2_2f [fitnumber]; T2_3 = T2_3f [fitnumber]; T2_4 = T2_4f [fitnumber]; A3_1 = A3_1f [fitnumber]; A3_2 = A3_2f [fitnumber]; A3_3 = A3_3f [fitnumber]; A3_4 = A3_4f [fitnumber]; T3_1 = T3_1f [fitnumber]; T3_2 = T3_2f [fitnumber]; T3_3 = T3_3f [fitnumber]; T3_4 = T3_4f [fitnumber]; D3_1 = D3_1f [fitnumber]; D3_2 = D3_2f [fitnumber]; D3_3 = D3_3f [fitnumber]; D3_4 = D3_4f [fitnumber]; t0 = t0f [fitnumber]; t1 = t1f [fitnumber]; t2 = t2f [fitnumber]; t3 = t3f [fitnumber]; t4 = t4f [fitnumber]; sprintf(filename, "%s/%s%i%s", datpath, datfile,fitnumber,".dat"); file = fopen("gnuplotfile", "w"); for (j=1; j<=10*t4; j++) fprintf (file, "%g %16.8g \n", (double)j/10., MODEL((double)j/10.) ); fclose(file); if ((fitnumber >= first)&&(fitnumber <= last)) { fprintf (pipe, "set term wxt noraise \n"); fprintf (pipe, "set title '<%s H=%.4f K=%.4f L=%.4f fcnred=%.2f>' \n", filename,Hf[fitnumber],Kf[fitnumber],Lf[fitnumber], FCNRED[fitnumber]); fprintf (pipe, "plot [:][:] '%s' using ($0):($2)/%g t '' w l lw 2 lt 1, " "'gnuplotfile' t '' w l lw 2 lt 2 \n", filename, MMON[fitnumber]); fflush(pipe); } } while (1==1); fclose (pipe); } else {/* makes 2D gnuplot files only if all data are involved: */ header(); file = fopen ("gnuplot2D_T1_1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A2_1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A2_1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A2_2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A2_2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A2_3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A2_3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A2_4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A2_4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T1_4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T1_4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T2_1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T2_1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T2_2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T2_2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T2_3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T2_3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_T2_4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], T2_4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S1*A1_1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S1f[tr(i)]*A1_1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S2*A1_2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S2f[tr(i)]*A1_2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S3*A1_3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S3f[tr(i)]*A1_3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S4*A1_4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S4f[tr(i)]*A1_4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_1*T1_1*S1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_1f[tr(i)]*T1_1f[tr(i)]*S1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_2*T1_2*S2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_2f[tr(i)]*T1_2f[tr(i)]*S2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_3*T1_3*S3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_3f[tr(i)]*T1_3f[tr(i)]*S3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_A1_4*T1_4*S4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], A1_4f[tr(i)]*T1_4f[tr(i)]*S4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S2-S1", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S2f[tr(i)] - S1f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S3-S2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S3f[tr(i)] - S2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S4-S3", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S4f[tr(i)] - S3f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S1-S4", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S1f[tr(i)] - S4f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file); file = fopen ("gnuplot2D_S4-S2", "w") ; fprintf (file, "%s", headerchar); for (i=first; i<=last; i++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(i)], Kf[tr(i)], S4f[tr(i)] - S2f[tr(i)]); if (i%31==0)fprintf (file, " \n");} fclose(file);} return;} /********************** function: read_fit_dat_files ***************/ /****************************** function: tr ***********************/ int tr (int i) {/*reorganising of the data file order: K will be monotonous*/ int j; if (i <= 310) j = i; else if (i <= 651) j = i + 310; else j = i - 341; return i; /* nem rendez at */ /*return j; atrendez */ } /****************************** function: tr ***********************/ /***************************** function: power *********************/ double power(double x, double p) {/* calculates x^p */ return exp(p * log(x));} /***************************** function: power *********************/ /********************** function: read_dat_files ***************/ void read_dat_files () { /* experimental counts will be read in: */ int ii, index=0, mtime, lineindex; char filename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; /* segedvaltozo */ int imon, idet; /* segedvaltozo egesztipusu */ char dummy[MAX_LEN]; FILE *file; for (ii=first; ii<=last; ii++){ sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); /* printf("%s \n", filename); */ printf("%i ", ii); fflush(stdout); file=fopen(filename,"r"); mtime=0; lineindex=0; MMONITOR = 0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) sscanf (line, "%s%f%f%f", dummy, &fH, &fK, &fL); else if (lineindex>12) {sscanf (line, "%i%i", &imon, &idet); mtime=mtime+1; index=index+1; Time[index]=mtime; H[index]=fH; K[index]=fK; L[index]=fL; mon[index]=imon; MMONITOR = MMONITOR + imon; det[index]=idet; if (maxdm<(double)idet/(double)imon) maxdm=(double)idet/(double)imon; /* printf("%i H= %f K=%f L=%f Time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],Time[index],mon[index],det[index]); */ } } MMON[ii] = MMONITOR/t4; fclose(file); } printf ("\nMAX(det/mon)= %f \n",maxdm); return;} /********************** function: read_dat_files ***************/ /**************** function: gnuplot_read_dat_files_t ***************/ void gnuplot_read_dat_files_t () { /* constructs 1D graphic files by gnuplot: file parameter: q; plot axis: time */ int ii, mtime, lineindex; char filename [MAX_LEN]; char gfilename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; char dummy[MAX_LEN]; FILE *file; FILE *pipe; /* H, K, L values will be read in: */ for (ii=first; ii<=last; ii++){ sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); /* here will be chosen: which graphics format will be used: */ /*P*/ /* sprintf(gfilename, "%s%03i%s", datfile,ii,".ps"); */ sprintf(gfilename, "%s%03i%s", datfile,ii,".png"); /* sprintf(gfilename, "%s%03i%s", datfile,ii,".gif"); */ printf("%i ", ii); fflush(stdout); file=fopen(filename,"r"); mtime=0; lineindex=0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) sscanf (line, "%s%f%f%f", dummy, &fH, &fK, &fL); /* printf("%i H= %f K=%f L=%f Time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],Time[index],mon[index],det[index]); */} fclose(file); /* gnuplot settings: */ pipe = popen("gnuplot -geometry 900x900+10+10", "w"); /* fprintf (pipe, "set logscale z \n"); */ fprintf (pipe, "set title '<%s H=%f K=%f L=%f>' \n", filename,fH,fK,fL); /* fprintf(pipe, "set terminal postscript color solid lw 2\n"); */ fprintf(pipe, "set terminal png \n"); /* fprintf(pipe, "set terminal gif \n"); */ fflush(pipe); fprintf(pipe, "set output '%s' \n", gfilename); fflush(pipe); /* gnuplot command: */ /* version constant scaling: */ /* fprintf (pipe, "plot [1:%f][0:%.0f] '%s' using ($0):($2)/'%f' t '' w l lw 2 lt 1 \n", t4, ceil(maxdm), filename, MMON [ii]); version auto scaling: */ fprintf (pipe, "plot [1:%f][0:] '%s' using ($0):($2)/'%f' t '' w l lw 2 lt 1 \n", t4, filename, MMON [ii]); fclose (pipe); } printf("\n"); return;} /**************** function: gnuplot_read_dat_files_t ***************/ /**************** function: gnuplot_read_dat_files_q ***************/ void gnuplot_read_dat_files_q () { /* constructs 2D graphic and gnuplot files by gnuplot: file parameter: time; plot axes: q */ int t, ii, lineindex, imon, idet; char filename [MAX_LEN]; char gfilename [MAX_LEN]; char line [MAX_LEN+1]; float fH, fK, fL; char dummy[MAX_LEN]; FILE *file; FILE *pipe; for (ii=first; ii<=last; ii++) { sprintf(filename, "%s/%s%i%s", datpath, datfile,ii,".dat"); file=fopen(filename,"r"); t=0; lineindex=0; if (file==NULL) terminate(" The desired file does not exist. "); while (fgets (line, MAX_LEN, file) != NULL) { lineindex=lineindex+1; if (lineindex==9) { sscanf (line, "%s%f%f%f", dummy, &fH, &fK, &fL); Hf[ii] = fH; Kf[ii] = fK; Lf[ii] = fL; /* printf("%i H= %f K=%f L=%f time=%i mon=%i det=%i \n ", index,H[index],K[index],L[index],time[index],mon[index],det[index]); */} else if (lineindex>12) { sscanf (line, "%i%i", &imon, &idet); t=t+1; COUNTSperMMON [ii][t] = idet/MMON [ii]; } } fclose(file);} for (t=1; t<=t4; t++){ sprintf(gfilename, "%s%04i%s", datfile,t,"t.gnu"); file = fopen (gfilename, "w"); for (ii=first; ii<=last; ii++){ fprintf (file, " %.4f \t %.4f \t %f \n", Hf[tr(ii)], Kf[tr(ii)], COUNTSperMMON [tr(ii)][t]); if (ii%31==0) fprintf (file, " \n");} fclose(file);} for (t=1; t<=t4; t++) { printf("%i ", t); fflush(stdout); /* here will be chosen: which graphics format will be used: */ /*P*/ /* sprintf(gfilename, "%s%04i%s", datfile,t,"t.ps"); */ sprintf( filename, "%s%04i%s", datfile,t,"t.png"); sprintf(gfilename, "%s%04i%s", datfile,t,"t.gnu"); /* sprintf(gfilename, "%s%04i%s", datfile,ii,"t.gif"); */ /* gnuplot settings: */ pipe = popen("gnuplot -geometry 900x900+10+10", "w"); /* fprintf (pipe, "set logscale zcb \n"); */ fprintf (pipe, "set title '<%s t=%i micro_s>' \n", filename, t); /* fprintf(pipe, "set terminal postscript color solid lw 2\n"); */ fprintf(pipe, "set terminal png \n"); /* fprintf(pipe, "set terminal gif \n"); */ fflush(pipe); fprintf(pipe, "set output '%s' \n", filename); fprintf(pipe, "set view 0,0 \n"); fflush(pipe); /* gnuplot command: */ /* version constant scaling: */ fprintf (pipe, "splot [][][0:%f] '%s' t '' with pm3d \n", ceil(maxdm), gfilename); /* version auto scaling: */ /* fprintf (pipe, "splot '%s' t '' with pm3d \n", gfilename); */ fclose (pipe); } printf("\n"); return;} /**************** function: gnuplot_read_dat_files_q ***************/ /********************** 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 and #include 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 (ERROR) terminate (" Floating point error detected.");} 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 ********************/