#define DATUM "2009_01_30" /******************************************X********** This is the file 'superfit.c.2009_01_30', genuine name: 'superfit.c'. The program 'superfit' 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. ************************************************** This program was originally written by Charlotte R. Anderson as 'fitspec.c' and 'fitoffspec.c'. The application of the minuit (CERN) was written by Andras Major as 'c-minuit'. The newest versions of the source code 'superfit.c', of the manual page 'superfit.1.gz', and of the c-minuit installation 'c-minuit' source file can freely be obtained from the web page of the Max Planck Institute for Metals Research (www.mf.mpg.de), Department Dosch, page "software". For more detail see the manual page 'superfit.1.gz'. ****************************************************/ /* run-time parameters (do not influence the results): */ #define SCREEN_WRITE 20 /* screen write out period in MINUIT fitting */ #define BORDER 30 /* substrate and vacuum width (%)in the nb plot */ /* end of run-time parameters */ /* for monitoring all errors: int ERROR = 0; for stopping at the first error: int ERROR = 1; */ int ERROR = 0; /* the next line may make the program slower but may be more accurate: */ #define __NO_MATH_INLINES #include #include #include #include #include #include #include /** comment out this line for DEC UNIX, please **/ #include #include "cfortran.h" #include "minuitfcn.h" #include "minuit.h" #define PI 3.1415926535897932384626433832795028841971693993751 #define PLANCK_H 6.62607554e-34 /* Planck's h in Js */ #define HBAR 1.0545726663e-34 /* in Js*/ #define M_N 1.67492861e-27 /* kg neutron mass */ #define MU_N 9.66237074e-27 /* neutron magnetic moment in J/T */ #define MAX_LEN 256 /* max. number of input characters */ #define NL 5000 /* max. number of layers */ #define NE 500000 /* max. number of data points */ #define GEN_PARAM 40 /* number of parameters in the header */ #define LAY_PARAM 7 /* number of parameters per layer */ #define MAXPIX 1500 /* number of pixels along the axis in plots */ #define EFFZERO 1e-40 /* my definition of zero */ #define IMAGZERO 1e-250 /* against underflow in complex division */ #define LIM_STOT 0.0001 /* error limit in defStot */ #define LIM_EXP 70. /* limit for the function exp */ #define LIM_defCn 0.5 /* see note at the top */ #define LIM_double 1.e-10 /* error of double numbers */ #define komplex complex double #define NFACTOR 1. /* multiplies all counts during read in */ /* definitions of structures that are needed in program: */ typedef struct {komplex x, y, z;} vector; /* 3D vector */ typedef struct {komplex m, n;} twoby1; /* 2D vector */ typedef struct {twoby1 Vm; twoby1 Vn;} svector; /*2D vector within 2D vector */ typedef struct {komplex a, b, c, d;} matrix; /* 2x2 matrix */ typedef struct {matrix Kx; matrix Ky; matrix Kz;} kmatrix; /* a set of three matrices */ typedef struct {matrix Sa; matrix Sb; matrix Sc; matrix Sd;} smatrix; /* 2x2 matrix with each component a 2x2 matrix = supermatrix */ /* used in main: */ char *name; /* generic name of the files involved */ char datfile [MAX_LEN]; /* name of file containing data to be fitted */ char minfile [MAX_LEN]; /* name of file defining minuit commands */ char minfilb [MAX_LEN]; /* backup of file defining minuit commands */ char ppmfile [MAX_LEN]; /* ppm 2D graphic file for the results */ char ppmfile_uu [MAX_LEN]; /* ppm 2D graphic file for the ++ results */ char ppmfile_ud [MAX_LEN]; /* ppm 2D graphic file for the +- results */ char ppmfile_du [MAX_LEN]; /* ppm 2D graphic file for the -+ results */ char ppmfile_dd [MAX_LEN]; /* ppm 2D graphic file for the -- results */ char ppmfile_au [MAX_LEN]; /* ppm 2D graphic file for the + asymmetry */ char ppmfile_ad [MAX_LEN]; /* ppm 2D graphic file for the - asymmetry */ char ascfile_uu [MAX_LEN]; /* ascii 2D file for the ++ results */ char ascfile_ud [MAX_LEN]; /* ascii 2D file for the +- results */ char ascfile_du [MAX_LEN]; /* ascii 2D file for the -+ results */ char ascfile_dd [MAX_LEN]; /* ascii 2D file for the -- results */ char psfile [MAX_LEN]; /* ps 1D file */ char pngfile [MAX_LEN]; /* png 1D file */ char nbpngfile [MAX_LEN]; /* png scattering length density profile file */ char nbpsfile [MAX_LEN]; /* png scattering length density profile file */ char simfile [MAX_LEN]; /* simulation file for the gnuplot */ char expfile [MAX_LEN]; /* experimental file for the gnuplot */ char modfile [MAX_LEN]; /* model file for the gnuplot */ char nbfile [MAX_LEN]; /* file for the gnuplot nb profile*/ /* experimental characteristics: */ double wavelen, bg_ampl, bg_alpha_i, bg_const, ampl; double P_analyse, P_polariz, N_analyse, N_polariz; vector Pvec, Povec; /* vectors defining polariser set-up */ matrix rho, rho0; /* analyser and polarizer density mat */ double blur_i, blur_f, blur_m; /* Gaussian blurrings in degrees */ double alpha_foot; /* angle limit for the footprint */ /* multilayer characteristics: */ int ndlayer; /* number of different layers in system */ int alllayer; /* total number of layers in system */ int npts, nrep[NL], pos[NL]; /* used to define proper order of m-l */ int translate[NL]; /* number of layer -> number of input layer */ char layername_in[NL][100]; /* reference names of the layers in minfile */ double B[NL], th[NL], fi[NL], d[NL], rough[NL]; double B_in[NL], th_in[NL], fi_in[NL], d_in[NL], rough_in[NL]; komplex nb[NL], nb_in[NL]; komplex nb_vacuum; double model_nob_1D[NE]; /* nb = scattering length density (m^-2) */ /* d = thickness of layer (m) */ /* rough = layer roughness (m) */ /* B = local magnetic field in m-l (T) */ /* th/fi = angles describing B orientation */ double chsix, chsiy, chsiz; /* correlation lengths */ double B_enhance; /* magnetic roughness enhancement factor (must be 1) */ double G_enhance; /* geometric roughness enhancement factor (must be 1) */ /* parameters used in calculations: */ double nb_M1; /* magnetic scattering length density for 1T */ double Fkonstant; /* constant in Fermi potential, see p. 18 of Squires */ double Hkonstant; /* constant for the calculation of the Hamiltonian */ double po[NL]; /* 0.5*Q-transfer perpendicular to m-l (z) */ matrix qmc2[NL], Rhat[NL]; smatrix Sm[NL]; /* the supermatrices divided by cosh(cimag())*/ matrix sm[NL]; /* to the supermatrices: cosh(cimag())*/ smatrix SmW[NL]; /* the supermatrices divided by cosh(cimag()) and multiplied by the Debye-Waller factor */ smatrix Stot; smatrix StotW; /* corrected by the Debye-Waller factors */ double detStot; /* cabs of the determinant of Stot */ komplex Cn[NL]; svector chicombi_in[NL], chicombi_out[NL]; /*wave function coefficients */ matrix Nn[NL]; smatrix Sm_in[NL]; matrix pm[NL]; matrix sm_in[NL]; int Parratt=0; /* 1 if Parratt algorithm will be used */ /* variables used in data file: */ double N[NE]; /* experimental counts */ double Nuu[NE], Nud[NE], Ndu[NE], Ndd[NE]; double BG[NE]; /* experimental background counts */ double alpha[NE]; /* experimental angles in degree: 1D */ double alpha_in[NE]; /* read-in experimental angles in degree: 1D */ double Err[NE]; /* experimental error */ double Erruu[NE], Errud[NE], Errdu[NE], Errdd[NE]; double monitor[NE]; /* monitor counts */ double monitoruu[NE], monitorud[NE], monitordu[NE], monitordd[NE]; double alpha_i[NE]; /* experimental grazing incidence angle in degree */ double alpha_f[NE]; /* experimental grazing final angle in degree */ double ai_offset; /* offset in the readdatfile function for the correction of the zero initial angle position */ double bg_factor; /* factor for manipulation of the experimental background: used in the readdatfile function */ int model_n; /* used in model_specular etc. */ double ai_delta, ai_min, ai_max; /* used in model_specular etc. */ double absorber[9]; /* absorbtion factors of the absorbers */ int agroup[NE]; /* group for absorber class */ /* used in fitting routine: */ double start[LAY_PARAM * NL], step[LAY_PARAM * NL]; double min[LAY_PARAM * NL], max[LAY_PARAM * NL]; char parname[LAY_PARAM * NL][MAX_LEN]; /* names of the MINUIT pars */ double fcn_end, fcnred_end; /* minfit results */ int n_par; /* minfit results */ int iexp = 0; /* number of experimental points */ int im = 1; /* number of fit parameters */ int imin, imax; /* fit limits */ double lower, higher; /* fit limits */ /* used in simulate: */ double mod [MAXPIX+1] [MAXPIX+1]; double out [MAXPIX+1] [MAXPIX+1], nor [MAXPIX+1] [MAXPIX+1]; /*in: blur*/ /* used for 2D plotting: */ int xnum, ynum; double sim[7][MAXPIX+1][MAXPIX+1]; /* linear scale */ double lsim[4][MAXPIX+1][MAXPIX+1]; /* log scale */ int si = 0; int fit_dim = 0; double function_min, function_max; float colour_min = 0.7, colour_max = 0.0; double x_min, x_max, y_min, y_max, degreeperxpixel, degreeperypixel; /******* prototype function titles ******************/ void readminfile (void); void readdatfile (void); void readpoldatfile (void); void rewrite (void); void simulate (void); void ppm (void); /* 2D plot */ void blur (double in [MAXPIX+1] [MAXPIX+1]); void blur1D (double in [MAXPIX+1]); /* NOT USED AT ALL */ void grafika (void); /* 1D plot */ void grafikaPNR (void); /*1D plot for PNR */ void nbplot (void); /* nb plot */ void colour (FILE *file, float h, float s, float v); double function (double x, double y); void getparameters (void); void test (void); void create_sample_files (void); void asymmetry (void); void log_2D (void); double untergrund (double alpha); double parratt(double fi); /* dwba model routines: */ double model_diffuse (double alpha_in_temp, double alpha_out_temp); double model_specular_basic (double alpha); double model_specular (double alpha); void model_specular_no_blurring (void); double footprint (double fi); double defOffSpec (double alpha_in_temp, double alpha_out_temp); komplex Sq_par (double alpha_in_temp, double alpha_out_temp); komplex s_nn_ (int n, int n_); void defCn (double alpha_in_temp, double alpha_out_temp); twoby1 deft0 (double effic, double gain); svector defchi0combi(twoby1 t0, twoby1 r0, double po_temp); void rearrange (void); matrix defNn (komplex nb_temp, double B_temp, double fi_temp, double th_temp); double defSpec (double alpha_in_temp); void defdmatrix (void); double defpo (double alpha_temp); void defSm (double po_temp); matrix defpm (komplex nb_temp, double B_temp, double fi_temp, double th_temp, double po_temp); matrix defqmc2 (komplex nb_temp, double B_temp, double fi_temp, double th_temp); matrix defHm (komplex nb_temp, double B_temp, double fi_temp, double th_temp); vector defBm (double B_temp, double fi_temp, double th_temp); vector defRbold (matrix Rhat_temp); matrix defRhat (double po_temp); smatrix defStot (void); smatrix defStotW (void); /* matrix manipulation routines: */ matrix expm(matrix M1); matrix cosm(matrix M1); matrix sinm(matrix M1); matrix rootm(matrix M1); matrix defeigenval(matrix M1); matrix invm(matrix M1); matrix malI0(komplex number); matrix matrixTIMEScomplex(matrix M1, komplex number); matrix squarem(matrix M1); matrix matrixTIMESmatrix(matrix M1, matrix M2); matrix matrixSUBmatrix(matrix M1, matrix M2); matrix matrixADDmatrix(matrix M1, matrix M2); matrix transposem(matrix M1); matrix conjm(matrix M1); komplex tracem(matrix M1); komplex detm(matrix M1); komplex det(komplex a, komplex b, komplex c, komplex d); /* vector manipulation routines: */ vector crossv(vector Vec1, vector Vec2); komplex dotv(vector Vec1, vector Vec2); vector addv(vector Vec1, vector Vec2); vector subv(vector Vec1, vector Vec2); vector conjv(vector Vec1); /* combinations of bras, kets, matrices, supermatrices and supervectors: */ svector matrixTIMESsvector (matrix M, svector S); komplex braket(twoby1 bra, twoby1 ket, matrix pot); twoby1 mtimesv(matrix M1, twoby1 ket); twoby1 add2by1(twoby1 t1, twoby1 t2); svector smatrixTIMESsvector(smatrix S, svector V); smatrix matrixTIMESsmatrix (matrix M, smatrix S); smatrix smatrixTIMESmatrix (smatrix S, matrix M); smatrix smatrixTIMEScomplex (smatrix S, komplex a); /* Miscellaneous functions: */ double to_radian (double x); double to_grad (double x); double alpha_to_qz (double alpha_i); double qz_to_alpha (double qz); komplex rcsin (komplex z); komplex rccos (komplex z); void terminate (char *text); void wait (char *text); void whistle (void); void message (char *text1, char *text2); void error (char *ref); int errorcount = 0; /* for the function error */ int mode; /* running mode */ int fcncount = 0; void apvector (char *text, vector x); void apsmatrix (char *text, smatrix x); void apmatrix (char *text, matrix x); void aptwoby1 (char *text, twoby1 x); int dwba_error = 0, dwba_ok = 0; int pixel; /* number of pixels in the result plots */ double k0; /* neutron wave number = 2 PI/lambda */ double harmonics = 1.; /* 1 for the beam and 2 for the first harmonics */ double w2harmonic = 0.; /* weight of the second harmonics 0 ... 1 */ int gnuplotversion; /* 0 if >= 4.0; >0 if <=3.7 */ int pivotoverflow = 1; /* counter for the function braket */ /* special fit features site#1: */ /* double Nb_m_s, NbTi_m_s; */ /******************************* main **************************/ int main (int argc, char *argv[]) { int idummy, maxcall; int dummy, i, j; char strategy[20], sdummy[MAX_LEN+1]; double ddummy; char *kontr; if (argc == 2){name=argv[1]; if (strcmp(name, "sample")==0) create_sample_files();} if (argc != 6) {message("\n You have made a mistake in the input,", "please try again using the format:\n"); message("", " superfit \n"); message (" = 0: test (values of the wave function will be printed):\n", " the data file will be read in if present but not used"); message (" = 1: 1D (specular) simulation: \n", " the data file will be read in if present but not used"); message (" = 11: 1D (specular) simulation using the Parratt algorithm:\n", " the data file will be read in if present but not used"); message (" = 2: 1D (specular) simulation: \n", " the data file will be read in and plotted"); message (" = 12: 1D (specular) simulation using the Parratt algorithm:\n", " the data file will be read in and plotted"); message (" = 3: 1D (specular) fit: \n", " the minuit command file will be rewritten"); message (" = 13: 1D (specular) fit using the Parratt algorithm: \n", " the minuit command file will be rewritten"); message (" = 22: 1D (specular) simulation: 4-segment PNR \n", " the data file will be read in and plotted"); message (" = 23: 1D (specular) fit: 4-segment PNR \n", " the minuit command file will be rewritten"); message (" = 5: 2D (specular and diffuse) simulation: \n", " the data file will be read in if present but not used"); message (" = 6: 2D (only the diffuse contribution) simulation: \n", " the data file will be read in if present but not used"); message (" = 7: 2D (specular and diffuse) fit: NOT READY YET\n", " the minuit command file will be rewritten"); message (" = 9:", "read-in test; the minuit command file will be rewritten\n"); message ("= the number of pixels in the plots\n", " 1D: 500; 2D: 100 (fast), 400 (quality).\n"); message (", : in degree units.","\n"); message ("For creating a command and a data sample file, please type:", "'superfit sample'.\n"); message ("Program version: superfit.c." DATUM, "" ); terminate("");} name=argv[1]; mode = strtod(argv[2], &kontr); if(kontr == argv[2]) terminate(" wrong second parameter\n"); pixel = floor(strtod(argv[3], &kontr)); if(kontr == argv[3]) terminate(" wrong third parameter \n"); lower = strtod(argv[4], &kontr); if(kontr == argv[4]) terminate(" wrong fourth parameter \n"); higher = strtod(argv[5], &kontr); if(kontr == argv[5]) terminate(" wrong fifth parameter \n"); y_min = lower; y_max = higher; degreeperypixel = (y_max-y_min)/(double)pixel; x_min = lower; x_max = higher; degreeperxpixel = (x_max-x_min)/(double)pixel; strcpy(datfile, name); strcat(datfile, ".dat"); strcpy(minfile, name); strcat(minfile, ".min"); strcpy(minfilb, name); strcat(minfilb, ".min~"); strcpy(ppmfile, name); strcat(ppmfile, ".ppm"); strcpy(ppmfile_uu, name); strcat(ppmfile_uu, "_uu.ppm"); strcpy(ppmfile_ud, name); strcat(ppmfile_ud, "_ud.ppm"); strcpy(ppmfile_du, name); strcat(ppmfile_du, "_du.ppm"); strcpy(ppmfile_dd, name); strcat(ppmfile_dd, "_dd.ppm"); strcpy(ascfile_uu, name); strcat(ascfile_uu, "_uu.asc"); strcpy(ascfile_ud, name); strcat(ascfile_ud, "_ud.asc"); strcpy(ascfile_du, name); strcat(ascfile_du, "_du.asc"); strcpy(ascfile_dd, name); strcat(ascfile_dd, "_dd.asc"); strcpy(ppmfile_au, name); strcat(ppmfile_au, "_au.ppm"); strcpy(ppmfile_ad, name); strcat(ppmfile_ad, "_ad.ppm"); strcpy(psfile, name); strcat(psfile, ".ps"); strcpy(pngfile, name); strcat(pngfile, ".png"); strcpy(nbpngfile, name); strcat(nbpngfile, "_nb.png"); strcpy(nbpsfile, name); strcat(nbpsfile, "_nb.ps"); strcpy(simfile, name); strcat(simfile, ".splot"); strcpy(expfile, name); strcat(expfile, ".eplot"); strcpy(modfile, name); strcat(modfile, ".mplot"); strcpy(nbfile, name); strcat(nbfile, ".nb"); message ("\n\n\n\n\n",""); message( " ______________________________________________________\n", "| |"); message (" | SUPERFIT |\n", "| |"); message (" | |\n", "| a fitting program |"); message (" | for specular and diffuse neutron reflectivity data |\n", "| obtained with polarized or non-polarized neutrons |"); message (" | on magnetic or non-magnetic samples |\n", "| |"); message (" | based on the combination of the Supermatrix Method |\n", "| and the Distorted Wave Born Approximation |"); message (" | |\n", "| |"); message (" | Max-Planck-Institut für Metallforschung |\n", "| Heisenbergstr. 3, D-70569 Stuttgart, Germany |"); message (" | |\n", "| (program version: " DATUM ") |"); message (" |______________________________________________________|", "\n\n"); printf ("The experimental file is called \"%s\". \n", datfile); printf ("The minuit command file is called \"%s\". \n", minfile); printf ("The minuit command file will be duplicated as \"%s\". \n", minfilb); printf ("The ++ ppm file will be called \"%s\". \n", ppmfile_uu); printf ("The +- ppm file will be called \"%s\". \n", ppmfile_ud); printf ("The -+ ppm file will be called \"%s\". \n", ppmfile_du); printf ("The -- ppm file will be called \"%s\". \n", ppmfile_dd); printf ("The ++ ascii file will be called \"%s\". \n", ascfile_uu); printf ("The +- ascii file will be called \"%s\". \n", ascfile_ud); printf ("The -+ ascii file will be called \"%s\". \n", ascfile_du); printf ("The -- ascii file will be called \"%s\". \n", ascfile_dd); printf ("The asymmetry+ ppm file will be called \"%s\". \n", ppmfile_au); printf ("The asymmetry- ppm file will be called \"%s\". \n", ppmfile_ad); printf ("The reference ppm file will be called \"%s\". \n", ppmfile); printf ("The postscript file will be called \"%s\". \n", psfile); printf ("The png file will be called \"%s\". \n", pngfile); printf ("The nb postscript file will be called \"%s\". \n", nbpsfile); printf ("The nb png file will be called \"%s\". \n", nbpngfile); printf ("The gnuplot simulation file will be called \"%s\". \n", simfile); printf ("The gnuplot experimental file will be called \"%s\". \n", expfile); printf ("The gnuplot model file will be called \"%s\". \n", modfile); printf ("The gnuplot nb file will be called \"%s\". \n", nbfile); printf ("\n"); /* first the necessary pre-calculations: */ nb_M1 = 2. * PI * MU_N * M_N / PLANCK_H / PLANCK_H; /* = magnetic scattering length density for 1 Tesla local magnetic field */ Fkonstant = 2. * PI * HBAR * HBAR / M_N; /* = constant in Fermi potential, see p. 18 of Squires */ if (mode == 9) printf ( "Magnetic scattering length density = %16.8e Tesla^-1 Angstrom^-2\n\n", nb_M1/1.e20); Hkonstant = 2. * M_N / (HBAR * HBAR); /* = constant for the calculation of the Hamiltonian */ /* first is ready */ /* second for 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");} /* second is ready*/ /* if (mode == 0) error ("main #1"); */ switch (mode){ case 0: /* test */ break; case 1: /* 1D simulation */ break; case 11: Parratt=1; /* 1D simulation */ break; case 2: /* 1D simulation + data plot */ break; case 12: Parratt=1; /* 1D simulation + data plot */ break; case 3: /* 1D fit */ fit_dim = 1; break; case 13: Parratt=1; /* 1D fit */ fit_dim = 1; break; case 22: /* 1D simulation of 4-segment PNR + data plot */ break; case 23: /* 1D fit of 4-segment PNR */ fit_dim = 1; break; case 5: /* 2D simulation */ break; case 6: /* 2D simulation (diffuse contributions only) */ break; case 7: /* 2D fit */ fit_dim = 2; break; case 9: /* rewrite .min */ break; default: terminate ("Illegal mode value.");} /* if (mode == 0) error ("main #2"); */ if (mode == 9) wait ("next: the minuit command file will be read in."); else if ((mode==0)||(mode==1)||(mode==11)||(mode==2)||(mode==12)||(mode==3)||(mode==13)||(mode==5)) message ("Next: the minuit command file will be read in.",""); readminfile(); if (mode == 9) wait ("next: the data file will be read in."); else if ((mode==1)||(mode==11)||(mode==2)||(mode==12)||(mode==3)||(mode==13)||(mode==5)) message ("Next: the data file will be read in.",""); if ((mode == 23)||(mode == 22)) readpoldatfile(); else readdatfile(); getparameters (); k0 = 2*PI/wavelen; if (alllayer +2 >= NL) terminate ("NL < number of layers !"); error ("begin"); if (mode == 9){ printf("The total number of layers is: %d\n", alllayer); printf("The number of different layers in the system is: %d\n", ndlayer); printf("The maximum number of fit parameters is: %d\n", im); wait("");} else message ("Please wait for the calculation.",""); /* if (mode == 0) error ("main #3"); */ switch (mode){ case 0: /* test */ test (); break; case 1: /* 1D simulation */ break; case 11: /*Parratt 1D simulation*/ break; case 2: /* 1D simulation + data plot */ break; case 12: /*Parratt 1D simulation + data plot*/ break; case 3: /* 1D fit */ break; case 13: /* Parratt 1D fit*/ break; case 22: /* 1D simulation 4-segment PNR + data plot*/ break; case 23: /* 1D fit 4-segment PNR */ break; case 5: /* 2D simulation */ ppm (); break; case 6: /* 2D simulation (diffuse contributions only) */ ppm (); break; case 7: /* 2D fit */ terminate ("Sorry: this mode (2D fit) is not implemented yet."); break; case 9: /* read-in test; rewrite .min */ rewrite (); break; default: terminate ("Illegal mode value.");} /* if (mode == 0) error ("main #4"); */ imin = 1; imax = iexp; for (i=1; i<=iexp; i++) if (lower > alpha[i]) imin = i+1; for (i=iexp; i>=1; i--) if (higher < alpha[i]) imax = i-1; if (imin < 1) imin = 1; if (imin > iexp) imin = iexp; if (imax > iexp) imax = iexp; if (imax < 1) imax = 1; printf ("Experimental points: number=%8d; fit range: %5d - %5d.\n", iexp, imin, imax); message("Please wait for the calculation.",""); if (fit_dim > 0) { MNINIT(5,6,7); /* do not change this line */ for (j=1; j<=im; j++) /* initializing parameters */ MNPARM(j, parname[j], start[j], step[j], min[j], max[j], dummy); maxcall = 4000; sprintf(strategy, "%s %i", "SIMPLEX", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ maxcall = 20000; sprintf(strategy, "%s %i", "MINIMIZE", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); /* fit */ /* maxcall = 10000; sprintf(strategy, "%s %i", "SEEK", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); */ /* fit */ /* maxcall = 8000; sprintf(strategy, "%s %i", "IMPROVE", maxcall); MNCOMD(minuitfcn, strategy, dummy, NULL); */ /* fit */ /* maxcall = 2000; sprintf(strategy, "%s %i", "MINIMIZE", maxcall); 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<=im; j++) /* fit results */ MNPOUT(j, sdummy, start[j], step[j], ddummy, ddummy, idummy); error ("ready"); rewrite ();} switch (mode){ case 0: /* test */ break; case 1: /* 1D simulation */ grafika (); break; case 11:/* Parratt 1D simulation */ grafika (); break; case 2: /* 1D simulation */ grafika (); break; case 12:/* Parratt 1D simulation */ grafika (); break; case 3: /* 1D fit */ grafika (); break; case 13:/* Parratt 1D fit */ grafika (); break; case 22:/*4-segment PNR 1D sim*/ grafikaPNR (); break; case 23:/*4-segment PNR 1D fit*/ grafikaPNR (); break; case 5: /* 2D simulation */ whistle (); break; case 6: /* 2D simulation (diffuse contributions only) */ whistle (); break; case 7: /* 2D fit */ break; case 9: /* rewrite .min */ break; default: terminate ("Illegal mode value.");} error("final"); if ((mode == 0)||(mode == 9)) #ifdef _FENV_H message ("The 'ISO C99 7.6: Floating-point environment '", "exists."); #else message ("The 'ISO C99 7.6: Floating-point environment '", "does not exist."); #endif if (fit_dim > 0) message("Please verify that no fit limit has been hit !", ""); return 0;} /****************************** main ***********************/ /***************************** fcn *****************************/ /* defines parameters' names and functions to be minimized */ void fcn(int npar,double* grad,double* fcnval,double* x,int iflag,void* futil){ int i, j; double model, dummy; wavelen = x[ 0]; P_polariz = x[ 1]; N_polariz = x[ 2]; P_analyse = x[ 3]; N_analyse = x[ 4]; ampl = x[ 5]; bg_const = x[ 6]; blur_i = x[ 7]; blur_f = x[ 8]; chsix = x[ 9]; chsiy = x[10]; chsiz = x[11]; B_enhance = x[12]; G_enhance = x[13]; alpha_foot = x[14]; ai_offset = x[15]; bg_factor = x[16]; bg_alpha_i = x[17]; bg_ampl = x[18]; w2harmonic = x[19]; absorber[1]= x[20]; absorber[2]= x[21]; absorber[3]= x[22]; absorber[4]= x[23]; absorber[5]= x[24]; /* special fit features site#2: */ /* Nb_m_s = x[25]; */ /* NbTi_m_s = x[26]; */ nb_vacuum = x[35] + I * x[36]; nb_in[0] = x[37] + I * x[38] - x[35]; nb[0] = nb_in[0] - x[35]; rough_in[0]= x[39]; rough[0] = x[39]; for(i=1; i<=ndlayer; i++){ nb_in[i] = x[(i*LAY_PARAM)+(GEN_PARAM -LAY_PARAM)] + x[(i*LAY_PARAM)+(GEN_PARAM+1-LAY_PARAM)] * I - x[35]; d_in [i] = x[(i*LAY_PARAM)+(GEN_PARAM+2-LAY_PARAM)]; rough_in[i] = x[(i*LAY_PARAM)+(GEN_PARAM+3-LAY_PARAM)]*d_in[i]; B_in[i] = x[(i*LAY_PARAM)+(GEN_PARAM+4-LAY_PARAM)]; if (x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)] > 180.) x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)] = x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)] - 360.; if (x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)] < -180.) x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)] = x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)] + 360.; fi_in[i] = x[(i*LAY_PARAM)+(GEN_PARAM+5-LAY_PARAM)]; if (x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)] > 180.) x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)] = x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)] - 360.; if (x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)] < -180.) x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)] = x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)] + 360.; th_in[i] = x[(i*LAY_PARAM)+(GEN_PARAM+6-LAY_PARAM)]; } nb_in[ndlayer+1] = 0.; /* uj */ blur_m = (blur_i + blur_f)/2.; /* special fit features site#3: */ /* d_in[1] = d_in[1] * Nb_m_s; */ /* d_in[2] = d_in[2] * NbTi_m_s; */ /* d_in[3] = d_in[3] * Nb_m_s; */ /* d_in[4] = d_in[4] * NbTi_m_s; */ /* d_in[5] = d_in[5] * Nb_m_s; */ for (i = 1; i <= iexp; i++) alpha[i] = alpha_in[i] - ai_offset; rearrange (); *fcnval=0; /* method of the least squares: for (j=imin; j<=imax; j++) {dummy = model(alpha[j])-N[j]; *fcnval+=dummy*dummy/(N[j]+1); }*/ /* maximum likelihood probability: */ if (mode == 23) { model_specular_no_blurring (); for (j=imin; j<=imax; j++) { dummy = Nuu[j]; model = model_specular(alpha[j]) * monitoruu[j] * absorber[agroup[j]]+BG[j]; if (model < 1.e-6) {message("Something is wrong: the model value is zero!","The"); message(datfile, "may be corrupted or in the "); message(minfile, "the background is set to zero strictly."); terminate("Please check it.");} *fcnval+=model-dummy; if (dummy>1.e-6) *fcnval+=dummy*log(dummy/model);} P_polariz = -P_polariz; model_specular_no_blurring (); for (j=imin; j<=imax; j++) { dummy = Ndu[j]; model = model_specular(alpha[j]) * monitordu[j] * absorber[agroup[j]]+BG[j]; if (model < 1.e-6) {message("Something is wrong: the model value is zero!","The"); message(datfile, "may be corrupted or in the "); message(minfile, "the background is set to zero strictly."); terminate("Please check it.");} *fcnval+=model-dummy; if (dummy>1.e-6) *fcnval+=dummy*log(dummy/model);} P_polariz = -P_polariz; P_analyse = -P_analyse; model_specular_no_blurring (); for (j=imin; j<=imax; j++) { dummy = Nud[j]; model = model_specular(alpha[j]) * monitorud[j] * absorber[agroup[j]]+BG[j]; if (model < 1.e-6) {message("Something is wrong: the model value is zero!","The"); message(datfile, "may be corrupted or in the "); message(minfile, "the background is set to zero strictly."); terminate("Please check it.");} *fcnval+=model-dummy; if (dummy>1.e-6) *fcnval+=dummy*log(dummy/model);} P_polariz = -P_polariz; model_specular_no_blurring (); for (j=imin; j<=imax; j++) { dummy = Ndd[j]; model = model_specular(alpha[j]) * monitordd[j] * absorber[agroup[j]]+BG[j]; if (model < 1.e-6) {message("Something is wrong: the model value is zero!","The"); message(datfile, "may be corrupted or in the "); message(minfile, "the background is set to zero strictly."); terminate("Please check it.");} *fcnval+=model-dummy; if (dummy>1.e-6) *fcnval+=dummy*log(dummy/model);} P_polariz = -P_polariz; P_analyse = -P_analyse; } else { if (fit_dim == 1) model_specular_no_blurring (); for (j=imin; j<=imax; j++) { dummy = N[j]; model = model_specular(alpha[j]) * monitor[j] * absorber[agroup[j]]+BG[j]; if (model < 1.e-6) {message("Something is wrong: the model value is zero!","The"); message(datfile, "may be corrupted or in the "); message(minfile, "the background is set to zero strictly."); terminate("Please check it.");} *fcnval+=model-dummy; if (dummy>1.e-6) *fcnval+=dummy*log(dummy/model); }} *fcnval = *fcnval*2.; /* print to screen: fit progress:*/ fcncount = fcncount + 1; if (fcncount%SCREEN_WRITE==0){ if (mode == 23) printf("fcn status: (%4i) fcn = %12.4f fcn_red = %12.4f \n", fcncount, *fcnval, *fcnval/(4.*(imax-imin+1)-npar)); else printf("fcn status: (%4i) fcn = %12.4f fcn_red = %12.4f \n", fcncount, *fcnval, *fcnval/((imax-imin+1)-npar));} if (iflag==3){fcn_end = *fcnval; fcnred_end = *fcnval/((imax-imin+1)-npar); if (mode == 23) fcnred_end = *fcnval/(4.*(imax-imin+1)-npar); n_par=npar;}} /**************************** fcn *********************************/ /************************** getparameters *****************************/ /* defines names of parameters in case of simulation only, not fit */ void getparameters(void) { int i; wavelen = start[ 1]; P_polariz = start[ 2]; N_polariz = start[ 3]; P_analyse = start[ 4]; N_analyse = start[ 5]; ampl = start[ 6]; bg_const = start[ 7]; blur_i = start[ 8]; blur_f = start[ 9]; chsix = start[10]; chsiy = start[11]; chsiz = start[12]; B_enhance = start[13]; G_enhance = start[14]; alpha_foot = start[15]; ai_offset = start[16]; bg_factor = start[17]; bg_alpha_i = start[18]; bg_ampl = start[19]; w2harmonic = start[20]; absorber[1] = start[21]; absorber[2] = start[22]; absorber[3] = start[23]; absorber[4] = start[24]; absorber[5] = start[25]; /* special fit features site#4: */ /* Nb_m_s = start[26]; */ /* NbTi_m_s = start[27]; */ nb_vacuum = start[36] + I * start[37]; nb_in[0] = start[38] + I * start[39] - start[36]; rough_in[0] = start[GEN_PARAM]; rough[0] = start[GEN_PARAM]; for(i=1; i<=ndlayer; i++){ nb_in[i] = start[(i-1)*LAY_PARAM+GEN_PARAM+1] + I * start[(i-1)*LAY_PARAM+GEN_PARAM+2] - start[36]; d_in [i] = start[(i-1)*LAY_PARAM+GEN_PARAM+3]; rough_in[i] = start[(i-1)*LAY_PARAM+GEN_PARAM+4]*d_in[i]; B_in[i] = start[(i-1)*LAY_PARAM+GEN_PARAM+5]; fi_in[i] = start[(i-1)*LAY_PARAM+GEN_PARAM+6]; th_in[i] = start[(i-1)*LAY_PARAM+GEN_PARAM+7];} nb_in[ndlayer+1] = 0.; /* uj */ blur_m = (blur_i + blur_f)/2.; /* special fit features site#5: */ /* d_in[1] = d_in[1] * Nb_m_s; */ /* d_in[2] = d_in[2] * NbTi_m_s; */ /* d_in[3] = d_in[3] * Nb_m_s; */ /* d_in[4] = d_in[4] * NbTi_m_s; */ /* d_in[5] = d_in[5] * Nb_m_s; */ for (i = 1; i <= iexp; i++) alpha[i] = alpha_in[i] - ai_offset; rearrange (); /* if (mode == 0) error ("getparameters"); */} /********************** getparameters *****************************/ /********************* rearrange *************************************/ /* rearranges the parameters into the actual layer format */ void rearrange (void){ int nlayer, i,j; int id = 0, il; nlayer = 0; nb[0]=nb_in[0]; /* substrat */ rough[0]=rough_in[0]; /* substrat */ /*printf( "\nrearrange: il= 0 nb=%8.2e rough=%8.2e (substrat) \n", nb[0], rough[0]);*/ while (id < ndlayer){ for (i = 1; i <= nrep[id+1]; i++) { for (j = 1; j <= pos[id+1]; j++){ il = nlayer + (i-1)*pos[id+1]+j; nb [il] = nb_in [id+j]; d [il] = d_in [id+j]; B [il] = B_in [id+j]; fi [il] = fi_in [id+j]; th [il] = th_in [id+j]; rough[il] = rough_in[id+j]; translate[il]= id + j; /*printf( "rearrange: il=%3d nb=%8.2e d=%8.2e rough=%8.2e B=% 8.2e\n", il, nb[il], d[il], rough[il], B[il]);*/ }} nlayer = nlayer + nrep[id+1]*pos[id+1]; /* number of layers (without substrat) */ id = id + pos[id+1]; } alllayer = nlayer; /* ez itt uj: */ translate[alllayer+1] = ndlayer + 1; nb[alllayer+1] = 0.; /* vacuum */ /*printf("rearrange: il=%3d nb=%8.2e (vacuum) \n\n", alllayer+1, nb[alllayer+1]); printf("rearrange: ndlayer=%4d nlayer=%4d \n", ndlayer, nlayer);*/ /* if (mode == 0) error ("rearrange"); */ } /********************** rearrange ************************/ /********************* model_specular ****************************/ double model_specular (double alpha) {/* alpha in degrees */ double dummy, coeff, mean, sum = 0., sumn = 0.; int i_mean, i, j; /* error("model_specular #0"); */ blur_m = (blur_i + blur_f)/2.; coeff = -4.0 * log(2.) * ai_delta * ai_delta / (blur_m * blur_m); /* that the Gaussian parameter = FWHM is */ mean = 1. + (alpha - ai_min)/ai_delta; i_mean = floor(mean); /* error("model_specular #1"); */ for (i=-2; i<=2; i++) {j = i + i_mean; if (j < 1) j = 1; if (j > model_n) j = model_n; dummy = exp(coeff * (i + i_mean - mean) * (i + i_mean - mean)); sum = sum + dummy * model_nob_1D[j]; sumn = sumn + dummy;} /* error("model_specular #2"); */ return (sum/sumn);} /********************* model_specular ****************************/ /********************* model_specular_no_blurring ********************/ void model_specular_no_blurring (void) { /* calculates the specular model for many points and put the values into model_nob_1D[i] */ double angle, firsthar, secondhar; int i; ai_min = alpha[1]; ai_max = alpha[iexp]; /*test line 13.07.2004*/ if (ai_min > lower) ai_min = lower; /*test line 13.07.2004*/ if (ai_max < higher) ai_max = higher; ai_delta = blur_m / 2.0; /* this 2.0 means that FWHM/2 is the density of the points */ if (ai_min <= ai_delta) ai_min = ai_delta; angle = ai_min - ai_delta; i = 0; while (angle <= ai_max) { i = i + 1; angle = angle + ai_delta; harmonics = 1.; firsthar = model_specular_basic(angle); harmonics = 2.; secondhar = model_specular_basic(angle); harmonics = 1.; model_nob_1D[i] = (1.-w2harmonic)*firsthar + w2harmonic*secondhar;} model_n = i; if (model_n > NE-1) {message("model_specular_no_blurring: the size of", "array is too small for this amount of data: please make NE larger."); terminate ("Solution for a more possible reason:" " please set the lower limit for the blur_i and blur_f to 0.01."); }} /********************* model_specular_no_blurring ********************/ /********************* model_specular_basic ****************************/ double model_specular_basic (double alpha) {/* alpha in degrees */ /* the Debye-Waller reduction of the specular reflectivity is included: only the roughness of the upmost layer is considered */ double dummy, DW_factor; dummy = rough[alllayer] * 2. * harmonics * k0 * sin(alpha * PI/180.); dummy = - dummy * dummy; DW_factor = dummy < -LIM_EXP ? 0 : exp(dummy); /* printf("footprint=%f defSpec=%f DW=%f\n", footprint(alpha),defSpec(alpha),DW_factor); */ if (Parratt==0) return (untergrund (alpha) + ampl * footprint (alpha) * defSpec (alpha) * DW_factor); else return (untergrund (alpha) + ampl * footprint (alpha) * parratt (alpha) * DW_factor); } /********************* model_specular_basic ****************************/ /************************* parratt ***************/ double parratt(double angle) { /* fi in Grad */ double ref, s2, ex, LAMBDA, LAMBDA2; int j, idummy; komplex k_ex, F[NL], f[NL], a[NL], R[NL], kdummy; LAMBDA = wavelen/harmonics; LAMBDA2 = LAMBDA * LAMBDA; s2=sin(angle * PI/180); s2 = s2 * s2; for (j=0; j<=alllayer+1; j++) /* original: f[j] = csqrt(s2-nb[j]*LAMBDA2/PI); modified for magnetic effects: */ f[j] = csqrt(s2 - (nb[j] + nb_M1 * B[j]) * LAMBDA2/PI); for (j=1; j<=alllayer+1; j++) { F[j] = (f[j]-f[j-1])/(f[j]+f[j-1]); /* Nevot-Croce roughness: */ ex=-8*PI*PI*rough[j]*rough[j]/LAMBDA2; k_ex=ex*f[j]*f[j-1]; idummy = cimag(k_ex)/(2.*PI); k_ex = k_ex - 2. * PI * idummy * I; k_ex = creal(k_ex) < -LIM_EXP ? 0 : cexp(k_ex); F[j]=k_ex*F[j];} for (j=1; j<=alllayer; j++){ kdummy = -I*PI/LAMBDA*f[j]*d[j]; idummy = cimag(kdummy)/(2.*PI); kdummy = kdummy -2.*PI*I*idummy; a[j] = cexp(kdummy);} /* kexp(kmul(kmul(k_i, dtok(-2*PI/lambda)), kmul(f[j], dtok(dl[j]))));} ----> version Stierle behind (16) */ /* kexp(kmul(kmul(k_i, dtok(-4*PI/lambda)), kmul(f[j], dtok(dl[j]))));} ----> version Parratt (7) */ a[alllayer+1] = 1.; R[0] = 0.; for (j=1; j<=alllayer+1; j++) {kdummy = 1.+ R[j-1] * F[j]; if (fabs(cimag(kdummy)) < IMAGZERO) kdummy = creal(kdummy); R[j] = a[j]*a[j]*a[j]*a[j]*(R[j-1]+ F[j])/kdummy;} ref = cabs(R[alllayer+1]); ref = ref * ref; return ref;} /************************* parratt ***************/ /********************* model_diffuse *******************************/ double model_diffuse(double alpha_in_temp, double alpha_out_temp){ /* alphas in degrees */ double theory, alpha_mean, speclim; alpha_mean = (alpha_in_temp + alpha_out_temp) / 2.; speclim = (higher - lower)/(double) pixel; /*** original version: if (mode == 3) if (fabs (alpha_in_temp - alpha_out_temp) < speclim) theory = untergrund (alpha_in_temp) + defOffSpec (alpha_in_temp, alpha_out_temp) + defSpec (alpha_mean); else theory = untergrund (alpha_in_temp) + footprint(alpha_in_temp) * defOffSpec (alpha_in_temp, alpha_out_temp); else theory = untergrund (alpha_in_temp) + footprint(alpha_in_temp) * defOffSpec (alpha_in_temp, alpha_out_temp); = end of the original version. actual version: ***/ if (mode != 6) if (fabs (alpha_in_temp - alpha_out_temp) < speclim) theory = untergrund (alpha_in_temp) + footprint(alpha_in_temp) * defOffSpec (alpha_in_temp, alpha_out_temp) + footprint(alpha_in_temp) * defSpec (alpha_mean); else theory = untergrund (alpha_in_temp) + footprint(alpha_in_temp) * defOffSpec (alpha_in_temp, alpha_out_temp); else theory = untergrund (alpha_in_temp) + footprint(alpha_in_temp) * defOffSpec (alpha_in_temp, alpha_out_temp); if (theory > 1.0) {theory = 0.; dwba_error++;} else dwba_ok++; return (theory);} /************************ model_diffuse ***********************************/ /********************** defOffSpec *************************/ /* defines off-specular reflectivity */ double defOffSpec(double alpha_in_temp, double alpha_out_temp) { /* alphas in degrees */ int i, j; komplex S_q, OffSpec, sum, factor; S_q = Sq_par(alpha_in_temp, alpha_out_temp); sum = 0; defCn(alpha_in_temp, alpha_out_temp); /* if (mode == 0) error ("defOffSpec:1"); */ /* perfect vertical correlation: for(i=0; i<=alllayer + 1; i++) for(j=0; j<=alllayer + 1; j++) sum = sum + S_q * Cn[i] * conj(Cn[j]); */ /* partial vertical correlation: */ for(i=0; i<=alllayer + 1; i++) for(j=0; j<=alllayer + 1; j++) { if (log10(cabs(1+Cn[i]))+log10(1+Cn[j]) > 299) {factor = 0.; terminate ("overflow in \"defOffSpec\" (too thick layers?) \n");} else factor = Cn[i] * conj(Cn[j]); sum = sum + s_nn_(i, j) * S_q * factor; /* sum = sum + s_nn_(i, j) * S_q * Cn[i] * conj(Cn[j]); */ /* printf("sum = %g + i*%g \n",creal(sum), cimag(sum)); */ /* printf("s_nn_(i, j) = %g + i*%g \n",creal(s_nn_(i,j)), cimag(s_nn_(i,j))); */ /* printf("sum = %g + i*%g \n",creal(sum), cimag(sum)); */ } /* vertically uncorrelated interfaces: for(j=0; j<=alllayer + 1; j++) sum = sum + S_q * Cn[j] * conj(Cn[j]); */ OffSpec = (k0*k0*k0/8./PI)*sum; return(OffSpec);} /************************* defOffSpec *********************************/ /************************* untergrund *********************************/ double untergrund (double alpha){/* calculates the background values for non-constant case */ /* alpha is in degrees */ /* the function is exponential decay */ double U; if (alpha > bg_alpha_i * LIM_EXP) U = bg_const; else U = bg_ampl * exp(-alpha/bg_alpha_i) + bg_const; return U;} /************************* untergrund *********************************/ /************************* s_nn_ ************************************/ komplex s_nn_ (int n, int n_) {/* correlation along z */ int i; double dd = 0, s, arg; if (n > n_) {i = n_; n_ = n; n = i;} for (i=n; i<=n_; i++) dd = dd + d[i]; /* s = exp(-0.5*dd*dd/(chsiz*chsiz)) / chsiz/sqrt(2.*PI); */ arg = -0.5*dd*dd/(chsiz*chsiz); s = arg < -LIM_EXP ? 0 : exp(arg); return s;} /************************* s_nn_ ************************************/ /************************* Sq_par ************************************/ /* correlations along the x and y directions */ komplex Sq_par(double alpha_in_temp, double alpha_out_temp){ /* alphas in degrees */ komplex K; double qx, arg; qx = k0*(cos(alpha_out_temp * PI/180.)-cos(alpha_in_temp * PI/180.)); arg = -0.5*qx*qx*chsix*chsix; arg = arg < -LIM_EXP ? 0 : exp(arg); K = sqrt(2.*PI)*chsix*arg; return (K);} /********************** Sq_par ***********************************/ /********************** defCn ******************************/ void defCn(double alpha_in_temp, double alpha_out_temp){ /* alphas in degrees */ int i, zero, zup, zdown; double po_in, po_out, konst, refm, refn; komplex konstant; matrix Rhat_in, Rhat_out, delVn; twoby1 r0_in, t0_in, r0_out, t0_out; svector chi0; chi0.Vm.m = 0.; chi0.Vm.n = 0.; chi0.Vn.m = 0.; chi0.Vn.n = 0.; po_in = defpo(alpha_in_temp); defSm (po_in); Rhat_in =defRhat(po_in); t0_in = deft0(P_polariz, N_polariz); r0_in = mtimesv(Rhat_in, t0_in); chicombi_in[alllayer+1] = defchi0combi (t0_in, r0_in, po_in); zero = 1; zup = 1; zdown = 1; refm = LIM_defCn * cabs(chicombi_in[alllayer+1].Vm.m); refn = LIM_defCn * cabs(chicombi_in[alllayer+1].Vm.n); for (i=alllayer; i>=0; i--) if (zero == 0) chicombi_in[i] = chi0; else {chicombi_in[i]= /* smatrixTIMESsvector (matrixTIMESsmatrix(sm[i],Sm[i]),chicombi_in[i+1]); */ /* matrixTIMESsvector(sm[i],smatrixTIMESsvector(SmW[i],chicombi_in[i+1])); */ matrixTIMESsvector(sm[i],smatrixTIMESsvector(Sm[i],chicombi_in[i+1])); if (fabs(detStot-1.) > 5) { /* message ("detStot /= 1","[pos #1]"); whistle(); */ if ((creal(chicombi_in[i].Vm.m) * creal(chicombi_in[i].Vn.m) > EFFZERO) || (cimag(chicombi_in[i].Vm.m) * cimag(chicombi_in[i].Vn.m) > EFFZERO)) {zup = 0; chicombi_in[i].Vm.m = 0.; chicombi_in[i].Vn.m = 0.;} if ((creal(chicombi_in[i].Vm.n) * creal(chicombi_in[i].Vn.n) > EFFZERO) || (cimag(chicombi_in[i].Vm.n) * cimag(chicombi_in[i].Vn.n) > EFFZERO)) {zdown = 0; chicombi_in[i].Vm.n = 0.; chicombi_in[i].Vn.n = 0.;} if (zup + zdown == 2) zero = 0; }} po_out = defpo(alpha_out_temp); defSm (po_out); Rhat_out =defRhat(po_out); t0_out = deft0(P_analyse, N_analyse); r0_out = mtimesv(Rhat_out, t0_out); chicombi_out[alllayer+1] = defchi0combi(t0_out, r0_out, po_out); zero = 1; zup = 1; zdown = 1; refm = LIM_defCn * cabs(chicombi_out[alllayer+1].Vm.m); refn = LIM_defCn * cabs(chicombi_out[alllayer+1].Vm.n); for(i=alllayer; i>=0; i--) if (zero == 0) chicombi_out[i] = chi0; else {chicombi_out[i]= /* smatrixTIMESsvector (matrixTIMESsmatrix(sm[i],Sm[i]),chicombi_out[i+1]); */ /* matrixTIMESsvector(sm[i],smatrixTIMESsvector(SmW[i],chicombi_out[i+1])); */ matrixTIMESsvector(sm[i],smatrixTIMESsvector(Sm[i],chicombi_out[i+1])); if (fabs(detStot-1.) > 1.5) { /* message ("detStot /= 1","[pos #2]" ); whistle(); */ if ((creal(chicombi_out[i].Vm.m) * creal(chicombi_out[i].Vn.m) > EFFZERO) || (cimag(chicombi_out[i].Vm.m) * cimag(chicombi_out[i].Vn.m) > EFFZERO)) {zup = 0; chicombi_out[i].Vm.m = 0.; chicombi_out[i].Vn.m = 0.;} if ((creal(chicombi_out[i].Vm.n) * creal(chicombi_out[i].Vn.n) > EFFZERO) || (cimag(chicombi_out[i].Vm.n) * cimag(chicombi_out[i].Vn.n) > EFFZERO)) {zdown = 0; chicombi_out[i].Vm.n = 0.; chicombi_out[i].Vn.n = 0.;} if (zup + zdown == 2) zero = 0; }} Nn[alllayer+1].a = 0; Nn[alllayer+1].b = 0; Nn[alllayer+1].c = 0; Nn[alllayer+1].d = 0; B[0] = 0.; fi[0] = 0.; th[0] = 0.; konstant = (4*PI*wavelen*wavelen)/(4*PI*PI); konst = sqrt (2/PI); for (i=alllayer; i >= 0; i--) { Nn[i] = defNn (nb[i], B_enhance * B[i], fi[i], th[i]); delVn = matrixTIMEScomplex(matrixSUBmatrix (Nn[i+1],Nn[i]),konstant); Cn[i] = konst * rough[i] * G_enhance * (braket(chicombi_in[i+1].Vm, chicombi_out[i+1].Vm, delVn)); /* if (i==0) { printf("chicombi_in[%d].Vm.m=%g+i*%g \n",i+1, creal(chicombi_in[i+1].Vm.m), cimag(chicombi_in[i+1].Vm.m)); printf("chicombi_in[%d].Vm.n=%g+i*%g \n",i+1, creal(chicombi_in[i+1].Vm.n), cimag(chicombi_in[i+1].Vm.n)); printf("chicombi_out[%d].Vm.m=%g+i*%g \n",i+1, creal(chicombi_out[i+1].Vm.m), cimag(chicombi_out[i+1].Vm.m)); printf("chicombi_out[%d].Vm.n=%g+i*%g \n",i+1, creal(chicombi_out[i+1].Vm.n), cimag(chicombi_out[i+1].Vm.n)); printf("rough[%d]=%g \n",i, rough[i]); printf("konst=%g \n",konst); printf("delVn.a=%g+i*%g \n", creal(delVn.a), cimag(delVn.a)); printf("delVn.b=%g+i*%g \n", creal(delVn.b), cimag(delVn.b)); printf("delVn.c=%g+i*%g \n", creal(delVn.c), cimag(delVn.c)); printf("delVn.d=%g+i*%g \n", creal(delVn.d), cimag(delVn.d)); printf("Cn[%d]=%g+i*%g \n",i, creal(Cn[i]), cimag(Cn[i])); printf("*Sm[%d].Sa.a=%g+i*%g \n",i, creal(Sm[i].Sa.a), cimag(Sm[i].Sa.a)); printf("*Sm[%d].Sb.a=%g+i*%g \n",i, creal(Sm[i].Sb.a), cimag(Sm[i].Sb.a)); printf("*Sm[%d].Sc.a=%g+i*%g \n",i, creal(Sm[i].Sc.a), cimag(Sm[i].Sc.a)); printf("*Sm[%d].Sd.a=%g+i*%g \n",i, creal(Sm[i].Sd.a), cimag(Sm[i].Sd.a)); printf("*Sm[%d].Sa.b=%g+i*%g \n",i, creal(Sm[i].Sa.b), cimag(Sm[i].Sa.b)); printf("*Sm[%d].Sb.b=%g+i*%g \n",i, creal(Sm[i].Sb.b), cimag(Sm[i].Sb.b)); printf("*Sm[%d].Sc.b=%g+i*%g \n",i, creal(Sm[i].Sc.b), cimag(Sm[i].Sc.b)); printf("*Sm[%d].Sd.b=%g+i*%g \n",i, creal(Sm[i].Sd.b), cimag(Sm[i].Sd.b)); printf("*sm[%d].a=%g+i*%g \n",i, creal(sm[i].a), cimag(sm[i].a)); }*/ }} /************************ defCn ***************************************/ /*************************** deft0 ****************************************/ /* defines the vector t0 using either the polariser of analyser efficiency*/ twoby1 deft0(double effic, double gain) { twoby1 t0; t0.m = sqrt(gain*(1 + effic)); t0.n = sqrt(gain*(1 - effic)); return(t0);} /************************ deft0 ***************************************/ /************************ defchi0combi *************************************/ /* defines the initial wave-function at the surface */ svector defchi0combi(twoby1 t0, twoby1 r0, double po_temp) { twoby1 chi0, chi0dash; svector chi0combi; chi0.m = t0.m + r0.m; chi0.n = t0.n + r0.n; chi0dash.m = (I*po_temp*t0.m)-(I*po_temp*r0.m); chi0dash.n = (I*po_temp*t0.n)-(I*po_temp*r0.n); chi0combi.Vm = chi0; chi0combi.Vn = chi0dash; if (mode == 0) { printf("defchi0combi: up t0=%+.2e%+.2ei r0=%+.2e%+.2ei\n", creal(t0.m), cimag(t0.m), creal(r0.m), cimag(r0.m)); printf("defchi0combi: down t0=%+.2e%+.2ei r0=%+.2e%+.2ei\n", creal(t0.n), cimag(t0.n), creal(r0.n), cimag(r0.n)); printf(" defchi0combi: Vm.m=%+.2e%+.2ei Vm.n=%+.2e%+.2ei \n", creal(chi0combi.Vm.m), cimag(chi0combi.Vm.m), creal(chi0combi.Vm.n), cimag(chi0combi.Vm.n)); printf(" defchi0combi: Vn.m=%+.2e%+.2ei Vn.n=%+.2e%+.2ei \n", creal(chi0combi.Vn.m), cimag(chi0combi.Vn.m), creal(chi0combi.Vn.n), cimag(chi0combi.Vn.n)); } return chi0combi;} /********************** defchi0combi ***********************************/ /*********************** defNn ***********************/ /* returns matrix defining nuclear and magnetic scattering density a la equation just under A5 of Nickel */ matrix defNn(komplex nb_temp, double B_temp, double fi_temp, double th_temp){ double nbmag; vector nmag_vec; matrix Nn_temp; nbmag = B_temp * nb_M1; nmag_vec=defBm(nbmag, fi_temp, th_temp); Nn_temp.a = nb_temp + nmag_vec.y; Nn_temp.b = nmag_vec.z - (I * nmag_vec.x); Nn_temp.c = nmag_vec.z + (I * nmag_vec.x); Nn_temp.d = nb_temp - nmag_vec.y; return Nn_temp;} /************************** end defNn ***************************/ /************************ defSpec *************************/ /* defines curly R from Adrian's paper */ double defSpec(double alpha_in_temp){ /* alpha in degrees */ komplex R_temp; double po_temp; matrix Rhat; double reflectivity; defdmatrix (); po_temp = defpo (alpha_in_temp); defSm (po_temp); Rhat = defRhat (po_temp); R_temp = tracem (matrixTIMESmatrix (rho, matrixTIMESmatrix (Rhat, matrixTIMESmatrix (rho0, transposem (conjm (Rhat)))))); reflectivity = cabs (R_temp); reflectivity = 4. * reflectivity; /* without this factor 4, for totalreflexion this function provides results, which are by a factor 4 too small -- reason: unknown */ return reflectivity;} /************************** defSpec ************************/ /*********************** defdmatrix ************************/ /* defines the density matrices from the polarization and analyser eff's given in input program. */ void defdmatrix(void){ Pvec.x = 0.; Pvec.y = P_analyse; Pvec.z = 0; Povec.x = 0.; Povec.y = P_polariz; Povec.z = 0; rho.a = (1. + Pvec.y) * N_analyse/2.; rho.b = (Pvec.z - I * Pvec.x) * N_analyse/2.; rho.c = (Pvec.z + I * Pvec.x) * N_analyse/2.; rho.d = (1. - Pvec.y) * N_analyse/2.; rho0.a = (1. + Povec.y) * N_polariz/2.; rho0.b = (Povec.z - I * Povec.x) * N_polariz/2.; rho0.c = (Povec.z + I * Povec.x) * N_polariz/2.; rho0.d = (1. - Povec.y) * N_polariz/2.;} /*********************** defdmatrix ********************************/ /************************** defpo ************************/ /* defines the component of the incoming wavevector that is perpendicular to the layers, po */ double defpo(double alpha_temp){ /* alpha in degrees */ return harmonics * k0 * sin(alpha_temp*PI/180.);} /*************************** defpo ***********************/ /************************ defRhat *********************/ /* now begins the definition of R-hat defined in equation 10 of Ru"hm et al. This will be used both in the definition of the specular part of the reflectivity and in the off specular part in order to define the wavefunction at the air-first layer interface, z = 0, equation 1 or Ru"hm et al. This wavefunction is needed so that the partial wavefunctions inside each layer may be calculated and used to define Cn of equation A6 in Nickel et al. */ matrix defRhat (double po_temp){ matrix curly, delta_1, Rhat_temp; matrix konst1, konst2 /* , konst1W, konst2W */; komplex ps, qsc; /* first calculate qsc, the critical wavenumber of the substrate */ Stot = defStot(); /* StotW = defStotW(); */ qsc=2*csqrt(PI*nb[0]); ps=csqrt((po_temp*po_temp)-(qsc*qsc)); /* the next two paragraphs are basically just defining delta and everything that is in the curly brackets of equation 10 in Ru"hm et al */ konst1 = matrixTIMEScomplex (Stot.Sc, (I / po_temp)); /* konst1W = matrixTIMEScomplex (StotW.Sc, (I / po_temp)); */ konst2 = matrixTIMEScomplex (Stot.Sb, (I * po_temp)); /* konst2W = matrixTIMEScomplex (StotW.Sb, (I * po_temp)); */ /* roughness is included: curly = matrixSUBmatrix (matrixTIMEScomplex (matrixSUBmatrix (StotW.Sd, konst1W), po_temp), matrixTIMEScomplex (matrixADDmatrix (StotW.Sa, konst2W), ps)); */ /* roughness is NOT included: */ curly = matrixSUBmatrix (matrixTIMEScomplex (matrixSUBmatrix (Stot.Sd, konst1), po_temp), matrixTIMEScomplex (matrixADDmatrix (Stot.Sa, konst2), ps)); delta_1= invm (matrixADDmatrix (matrixTIMEScomplex (matrixADDmatrix (Stot.Sd, konst1), po_temp), matrixTIMEScomplex (matrixSUBmatrix (Stot.Sa, konst2), ps))); Rhat_temp = matrixTIMESmatrix (delta_1, curly); return Rhat_temp;} /************************** defRhat *************************/ /********************** defStot ***********************/ /*defines the total supermatrix for the whole multilayer system*/ smatrix defStot() { int i; smatrix Stemp, Sttemp; /* double control; */ Stemp = Sm[alllayer]; for(i=alllayer-1; i>=1; i--){ /* ez maradt valtozatlanul */ Sttemp.Sa = matrixADDmatrix (matrixTIMESmatrix (Sm[i].Sa, Stemp.Sa), matrixTIMESmatrix (Sm[i].Sb, Stemp.Sc)); Sttemp.Sb = matrixADDmatrix (matrixTIMESmatrix (Sm[i].Sa, Stemp.Sb), matrixTIMESmatrix (Sm[i].Sb, Stemp.Sd)); Sttemp.Sc = matrixADDmatrix (matrixTIMESmatrix (Sm[i].Sc, Stemp.Sa), matrixTIMESmatrix (Sm[i].Sd, Stemp.Sc)); Sttemp.Sd = matrixADDmatrix (matrixTIMESmatrix (Sm[i].Sc, Stemp.Sb), matrixTIMESmatrix (Sm[i].Sd, Stemp.Sd)); Stemp = Sttemp; /* apsmatrix("defStot: Sm[i] = ",Sm[i]); apsmatrix("defStot: Stemp = ",Stemp); if (mode == 0){ control = cabs(detm(matrixSUBmatrix(matrixTIMESmatrix(Stemp.Sa,Stemp.Sd), matrixTIMESmatrix(Stemp.Sb,Stemp.Sc))));}*/ } /* if (mode == 0) error ("defStot"); */ detStot = cabs(detm(matrixSUBmatrix(matrixTIMESmatrix(Stemp.Sa,Stemp.Sd), matrixTIMESmatrix(Stemp.Sb,Stemp.Sc)))); return Stemp;} /************************** defStot *************************/ /********************** defStotW ***********************/ /*defines the total supermatrix for the whole multilayer system corrected by the Debye-Waller factors */ smatrix defStotW() { int i; smatrix Stemp, Sttemp; Stemp = SmW[alllayer]; for(i=alllayer-1; i>=1; i--){ /* ez maradt valtozatlanul */ Sttemp.Sa = matrixADDmatrix (matrixTIMESmatrix (SmW[i].Sa, Stemp.Sa), matrixTIMESmatrix (SmW[i].Sb, Stemp.Sc)); Sttemp.Sb = matrixADDmatrix (matrixTIMESmatrix (SmW[i].Sa, Stemp.Sb), matrixTIMESmatrix (SmW[i].Sb, Stemp.Sd)); Sttemp.Sc = matrixADDmatrix (matrixTIMESmatrix (SmW[i].Sc, Stemp.Sa), matrixTIMESmatrix (SmW[i].Sd, Stemp.Sc)); Sttemp.Sd = matrixADDmatrix (matrixTIMESmatrix (SmW[i].Sc, Stemp.Sb), matrixTIMESmatrix (SmW[i].Sd, Stemp.Sd)); Stemp = Sttemp;} /* if (mode == 0) error ("defStotW"); */ return Stemp;} /************************** defStotW *************************/ /*********************** defSm *************************/ /* defines supermatrix for layer with entered d, nb, B, fi & th and for incoming angle defined by po_temp The supermatrices that are defined here are kept in the global variable, Sm[layer], and combinations of the matrices are used both in defStot, for defining Rhat and also in defOffSpec for */ void defSm (double po_temp) { matrix pmd, cos_pmd, sin_pmd, lambda, pm1, pm2; int i; matrix contr, cosh_pmdi; double aa,ba,thba,thbd,sinaa,sinad,cosaa,cosad,coshba,coshbd,sinhba,sinhbd; double ad,bd; komplex rcsina, rcsind, rccosa, rccosd; /* pm[ndlayer+1] = defpm(nb_in[ndlayer+1], B_in[ndlayer+1], fi_in[ndlayer+1], th_in[ndlayer+1], po_temp); */ pm[ndlayer+1] = defpm(nb_vacuum, 0, 0, 0, po_temp); for (i=1; i<=ndlayer; i++){/* loop over the different layers */ pm[i] = defpm(nb_in[i], B_in[i], fi_in[i], th_in[i], po_temp); pmd = matrixTIMEScomplex(pm[i], d_in[i]); /* error("defSm #1"); */ /* the sin and cos of pmd will be calculated: */ lambda = defeigenval(pmd); /* error("defSm #2"); */ if (mode == 0) { printf ("pmd.a= %+.4e%+.4ei pmd.b= %+.4e %+.4ei \n", creal(pmd.a), cimag(pmd.a), creal(pmd.b), cimag(pmd.b)); printf ("pmd.c= %+.4e%+.4ei pmd.d= %+.4e %+.4ei \n", creal(pmd.c), cimag(pmd.c), creal(pmd.d), cimag(pmd.d)); printf ("lambda.a= %+.4e%+.4ei lambda.d= %+.4e %+.4ei \n", creal(lambda.a), cimag(lambda.a), creal(lambda.d), cimag(lambda.d));} if (cabs (lambda.a-lambda.d) < LIM_STOT){ aa = creal(lambda.a); ba = cimag(lambda.a); bd = -8.88e-88; thba = tanh(ba); thbd = -8.88e-88; sinaa = sin(aa); sinad = -8.88e-88; cosaa = cos(aa); cosad = -8.88e-88; rcsina = sinaa + cosaa * thba * I; rcsind = -8.88e-88; rccosa = cosaa - sinaa * thba * I; rccosd = -8.88e-88; /* if overflow: */ if (ba < -LIM_EXP) terminate ("defSm: argument of cosh < -LIM_EXP (1)"); if (ba > LIM_EXP) {coshba = 0.; sinhba = 0.;} else {coshba = cosh(ba); sinhba =sinh(ba);} coshbd = -8.88e-88; sinhbd = -8.88e-88; cos_pmd = matrixADDmatrix(malI0(rccosa+lambda.a*rcsina), matrixTIMEScomplex(pmd,-rcsina)); sin_pmd = matrixADDmatrix(malI0(rcsina-lambda.a*rccosa), matrixTIMEScomplex(pmd, rccosa)); cosh_pmdi = matrixADDmatrix(malI0(coshba-lambda.a*sinhba), matrixTIMEScomplex(pmd, sinhba));} else { aa = creal(lambda.a);ad = creal(lambda.d); ba = cimag(lambda.a);bd = cimag(lambda.d); thba = tanh(ba); thbd = tanh(bd); sinaa = sin(aa); sinad = sin(ad); cosaa = cos(aa); cosad = cos(ad); rcsina = (sinaa + I*cosaa * thba); rcsind = (sinad + I*cosad * thbd); rccosa = (cosaa - I*sinaa * thba); rccosd = (cosad - I*sinad * thbd); /* if overflow: */ if (ba < -LIM_EXP) terminate ("defSm: argument of cosh < -LIM_EXP (2)"); if (bd < -LIM_EXP) terminate ("defSm: argument of cosh < -LIM_EXP (3)"); if (ba > LIM_EXP) coshba = 0.; else coshba = cosh(ba); if (bd > LIM_EXP) coshbd = 0.; else coshbd = cosh(bd); cos_pmd= matrixADDmatrix(malI0((rccosa*lambda.d-lambda.a*rccosd)/ (lambda.d-lambda.a)), matrixTIMEScomplex(pmd,((rccosd-rccosa)/(lambda.d-lambda.a)))); sin_pmd=matrixADDmatrix(malI0((rcsina*lambda.d-lambda.a*rcsind)/ (lambda.d-lambda.a)), matrixTIMEScomplex(pmd,((rcsind-rcsina)/(lambda.d-lambda.a)))); cosh_pmdi=matrixADDmatrix(malI0((coshba*lambda.d-lambda.a*coshbd)/ (lambda.d-lambda.a)), matrixTIMEScomplex(pmd,((coshbd-coshba)/(lambda.d-lambda.a))));} /* the sin and cos of pmd was calculated */ if (mode == 0) { printf("control: thba=%.2e thbd=%.2e coshba=%.2e coshbd=%.2e\n", thba, thbd, coshba, coshbd); printf("control: sinaa=%.2e cosaa=%.2e sinad=%.2e cosad=%.2e\n", sinaa, cosaa, sinad, cosad); printf("control: rcsina=%.2e+%.2ei rccosa=%.2e+%.2ei\n", creal(rcsina), cimag(rcsina), creal(rccosa), cimag(rccosa)); contr = matrixADDmatrix(squarem(sin_pmd),squarem(cos_pmd)); printf("control_pmd: %+.4e%+.4ei %+.4e%+.4ei\n " " %+.4e%+.4ei %+.4e%+.4ei\n", creal(contr.a), cimag(contr.a), creal(contr.b), cimag(contr.b), creal(contr.b), cimag(contr.b), creal(contr.d), cimag(contr.d) ); contr = sin_pmd; printf("sin_pmd: %+.4e%+.4ei %+.4e%+.4ei\n " " %+.4e%+.4ei %+.4e%+.4ei\n", creal(contr.a), cimag(contr.a), creal(contr.b), cimag(contr.b), creal(contr.b), cimag(contr.b), creal(contr.d), cimag(contr.d) ); contr = cos_pmd; printf("cos_pmd: %+.4e%+.4ei %+.4e%+.4ei\n " " %+.4e%+.4ei %+.4e%+.4ei\n", creal(contr.a), cimag(contr.a), creal(contr.b), cimag(contr.b), creal(contr.b), cimag(contr.b), creal(contr.d), cimag(contr.d) ); contr = cosh_pmdi; printf("cosh_pmdi: %+.4e%+.4ei %+.4e%+.4ei\n " " %+.4e%+.4ei %+.4e%+.4ei\n", creal(contr.a), cimag(contr.a), creal(contr.b), cimag(contr.b), creal(contr.b), cimag(contr.b), creal(contr.d), cimag(contr.d) );} Sm_in[i].Sa=cos_pmd; if (cabs(detm(pm[i]))<1.e-20) terminate ("defSm: cabs(detm(pm[i]))<1.e-20)"); else { Sm_in[i].Sb=matrixTIMESmatrix(matrixTIMEScomplex(invm(pm[i]), +1.), sin_pmd); Sm_in[i].Sc=matrixTIMESmatrix(matrixTIMEScomplex( pm[i], -1.), sin_pmd); } Sm_in[i].Sd=cos_pmd; sm_in[i] = cosh_pmdi;} /* ez is maradt valtozatlanul */ for (i=1; i<=alllayer; i++) {/* translates to the real layer system and corrects for the Debye-Waller exponential: ONLY the real parts of the pm normal momenta are included */ pm1.a = creal(pm[translate[i]].a); pm1.b = creal(pm[translate[i]].b); pm1.c = creal(pm[translate[i]].c); pm1.d = creal(pm[translate[i]].d); pm2.a = creal(pm[translate[i+1]].a); pm2.b = creal(pm[translate[i+1]].b); pm2.c = creal(pm[translate[i+1]].c); pm2.d = creal(pm[translate[i+1]].d); Sm[i] = Sm_in[translate[i]]; sm[i] = sm_in[translate[i]];} /* if (mode == 0) error ("defSm"); */ } /************************** defSm **** ***********************/ /*********************** defpm ***********************/ /*returns matrix pm (see page 16075 or Ru"hm et al)*/ matrix defpm (komplex nb_temp, double B_temp, double fi_temp, double th_temp, double po_temp) { matrix pm, qmc2_temp; qmc2_temp=defqmc2(nb_temp, B_temp, fi_temp, th_temp); pm=rootm(matrixSUBmatrix(malI0(po_temp*po_temp),qmc2_temp)); return pm;} /************************** defpm *************************/ /*********************** defqmc2 **********************/ /* returns matrix q(mc)^2, the eigenvalues of which are the critical wave numbers of reflections for neutrons with spin parallel and antiparallel to Bm */ matrix defqmc2(komplex nb_temp, double B_temp, double fi_temp, double th_temp){ matrix qmc2_temp, Hm; Hm=defHm(nb_temp, B_temp, fi_temp, th_temp); /* define Hamiltonian*/ qmc2_temp.a=Hkonstant*Hm.a; qmc2_temp.b=Hkonstant*Hm.b; qmc2_temp.c=Hkonstant*Hm.c; qmc2_temp.d=Hkonstant*Hm.d; return qmc2_temp;} /************************** defqmc2 ***************************/ /*********************** defHm ***********************/ /* returns matrix defining Hamiltonian within each layer equal for identical layers */ matrix defHm(komplex nb_temp, double B_temp, double fi_temp, double th_temp){ komplex VN_temp; vector Bm; matrix Hm; /* = constant in Fermi potential, see p.18 of Squires */ VN_temp = Fkonstant * nb_temp; Bm = defBm (B_temp, fi_temp, th_temp); Hm.a = VN_temp - MU_N * Bm.y; Hm.b = -1.* MU_N * (Bm.z-I*Bm.x); Hm.c = -1.* MU_N * (Bm.z+I*Bm.x); Hm.d = VN_temp + MU_N * Bm.y; return Hm;} /************************** defHm ***************************/ /*********************** defBm ***********************/ /* defines Bx, By and Bz, which are constant for each layer set */ vector defBm (double B_temp, double fi_temp, double th_temp) { vector Bm; double th_temp_rad, fi_temp_rad; th_temp_rad = th_temp * PI / 180.; fi_temp_rad = fi_temp * PI / 180.; Bm.x = B_temp * cos(th_temp_rad) * sin(fi_temp_rad); Bm.y = B_temp * cos(th_temp_rad) * cos(fi_temp_rad); Bm.z = B_temp * sin(th_temp_rad); return Bm;} /************************** defBm *********************************/ /************************ defRbold ******************************/ /* returns a vector = Trace (Rhat times Pauli matrices) */ vector defRbold(matrix Rhat_temp) { kmatrix sigma, Rsigma; vector Rbold_temp; sigma.Kz.a = 0.; sigma.Kz.b = 1.; sigma.Kz.c = 1.; sigma.Kz.d = 0.; sigma.Kx.a = 0.; sigma.Kx.b = -1.*I; sigma.Kx.c = I; sigma.Kx.d = 0.; sigma.Ky.a = 1.; sigma.Ky.b = 0.; sigma.Ky.c = 0.; sigma.Ky.d = -1.; Rsigma.Kx = matrixTIMESmatrix (Rhat_temp, sigma.Kx); Rsigma.Ky = matrixTIMESmatrix (Rhat_temp, sigma.Ky); Rsigma.Kz = matrixTIMESmatrix (Rhat_temp, sigma.Kz); Rbold_temp.x = tracem (Rsigma.Kx); Rbold_temp.y = tracem (Rsigma.Ky); Rbold_temp.z = tracem (Rsigma.Kz); return Rbold_temp;} /********************* defRbold ****************************/ /*********************** foot print ***************/ double footprint (double fi){ /* fi, alpha_foot in Grad */ double bel; /* the simple approach: if (fi < alpha_foot) bel = sin(fi * PI/180.) / sin(alpha_foot * PI/180.); else bel = 1; bel = (1+bel)/2; */ bel = sin(fi*PI/180.) / sin(alpha_foot*PI/180.) * sqrt(log(2.)); bel = bel > LIM_EXP ? 1 : erf (bel); return bel;} /*********************** foot print ***************/ /* matrix manipulations: */ /****************************** expm ***************************/ matrix expm (matrix M1){ komplex fa, fb; matrix M3, lambda; double n; n = cabs(M1.a) + cabs(M1.b) + cabs(M1.c) + cabs(M1.d); /* printf ("expm: fabs = %.4e \n", n); */ if (fabs(n) < EFFZERO) { M3 = matrixADDmatrix( malI0(1.), M1);} else {lambda = defeigenval(M1); if (cabs (lambda.a-lambda.d) < EFFZERO){ if (creal(lambda.a)<-LIM_EXP) fa = 0.; else fa = cexp(lambda.a); fb = fa; M3 = matrixADDmatrix(malI0(fa-lambda.a*fb), matrixTIMEScomplex(M1,fb));} else { if (creal(lambda.a)<-LIM_EXP) fa = 0.; else fa = cexp(lambda.a); if (creal(lambda.d)<-LIM_EXP) fb = 0.; else fb = cexp(lambda.d); M3 = matrixADDmatrix(malI0((fa*lambda.d-lambda.a*fb)/(lambda.d-lambda.a)), matrixTIMEScomplex(M1,((fb-fa)/(lambda.d-lambda.a))));}} return M3;} /***************************** expm *************************/ /****************************** cosm ***************************/ matrix cosm (matrix M1){ komplex fa, fb; matrix M3, lambda; lambda = defeigenval(M1); if (cabs (lambda.a-lambda.d) < EFFZERO){ /* fa = rccos(lambda.a); fb = -rcsin(lambda.a); */ fa = ccos(lambda.a); fb = -csin(lambda.a); M3 = matrixADDmatrix(malI0(fa-lambda.a*fb), matrixTIMEScomplex(M1,fb));} else { /* fa = rccos(lambda.a); fb = rccos(lambda.d); */ fa = ccos(lambda.a); fb = ccos(lambda.d); M3 = matrixADDmatrix(malI0((fa*lambda.d-lambda.a*fb)/(lambda.d-lambda.a)), matrixTIMEScomplex(M1,((fb-fa)/(lambda.d-lambda.a))));} return M3;} /***************************** cosm *************************/ /****************************** sinm ***************************/ matrix sinm (matrix M1){ komplex fa, fb; matrix M3, lambda; lambda = defeigenval(M1); if (cabs (lambda.a-lambda.d) < EFFZERO){ /* fa = rcsin(lambda.a); fb = rccos(lambda.a); */ fa = csin(lambda.a); fb = ccos(lambda.a); M3 = matrixADDmatrix(malI0(fa-lambda.a*fb), matrixTIMEScomplex(M1,fb));} else { /* fa = rcsin(lambda.a); fb = rcsin(lambda.d); */ fa = csin(lambda.a); fb = csin(lambda.d); M3 = matrixADDmatrix(malI0((fa*lambda.d-lambda.a*fb)/(lambda.d-lambda.a)), matrixTIMEScomplex(M1,((fb-fa)/(lambda.d-lambda.a))));} return M3;} /***************************** sinm *************************/ /****************************** rootm ***************************/ matrix rootm(matrix M1){/*returns matrix which is root of matrix M1*/ matrix lambda, M2, M3; komplex x, y, z, w; double n; /*matrix matrixx;*/ n = cabs(M1.a) + cabs(M1.b) + cabs(M1.c) + cabs(M1.d); if (fabs(n) < EFFZERO) terminate ("rootm: a matrix szingularis !"); M2.a = M1.a / n; M2.b = M1.b / n; M2.c = M1.c / n; M2.d = M1.d / n; n = sqrt (n); lambda=defeigenval(M2); x = csqrt(lambda.a); y = csqrt(lambda.d); z = x + y; w = x * y; M3.a = n * (M2.a + w) / z; M3.b = n * M2.b / z; M3.c = n * M2.c / z; M3.d = n * (M2.d + w) / z; /*matrixx = matrixSUBmatrix(squarem(M3),M1); if (cabs(matrixx.a)>.0099){ printf ("matrix: %12.2e %12.2e %12.2e %12.2e \n", cabs(M1.a), cabs(M1.b), cabs(M1.c), cabs(M1.d)); printf ("l_1, l_2, z, w: %12.2e %12.2e %12.2e %12.2e \n", cabs(lambda.a), cabs(lambda.d), cabs(z), cabs(w)); printf ("rootm: %12.2e %12.2e %12.2e %12.2e \n", cabs(M3.a), cabs(M3.b), cabs(M3.c), cabs(M3.d)); printf ("diff: %12.2e %12.2e %12.2e %12.2e \n", cabs(matrixx.a), cabs(matrixx.b), cabs(matrixx.c), cabs(matrixx.d));}*/ return M3; } /***************************** rootm *************************/ /****************************** defeigenval *************************/ /* defines eigenvalues of matrix M1 in form of matrix, (1,1)=lambda1, (1,2)=0, (2,1)=0, (2,2)=lambda2. */ matrix defeigenval (matrix M1){ matrix lambda, M2; komplex TrM, DtM, dummy; double n; n = cabs(M1.a) + cabs(M1.b) + cabs(M1.c) + cabs(M1.d); if (fabs(n) < EFFZERO) terminate ("defeigenval: the matrix is a nullmatrix:\n" "I cannot calculate the eigenvalues !"); M2.a = M1.a / n; M2.b = M1.b / n; M2.c = M1.c / n; M2.d = M1.d / n; DtM=detm(M2); /* if(cabs(DtM) max) max = cabs(b); if (cabs(c) > max) max = cabs(c); if (cabs(d) > max) max = cabs(d); if (max > 1.) {aa = a/max; bb = b/max; cc = c/max; dd = d/max;} else {aa = a; bb = b; cc = c; dd = d; max = 1.;} return max * max * ((aa * dd) - (bb * cc));} /************************** det ********************************/ /************************ transposem ******************/ matrix transposem (matrix M1){/*returns the transpose matrix of matrix M1*/ matrix M3; M3.a = M1.a; M3.b = M1.c; M3.c = M1.b; M3.d = M1.d; return M3;} /************************** transposem ************************/ /************************ conjm ******************/ matrix conjm (matrix M1){/*returns the conjugate matrix of matrix M1*/ matrix M3; M3.a = conj (M1.a); M3.b = conj (M1.b); M3.c = conj (M1.c); M3.d = conj (M1.d); return M3;} /************************** conjm ************************/ /* Functions for vector manipulation: */ /************************ crossv ******************/ /* Vec1 x Vec2 */ vector crossv(vector Vec1, vector Vec2){ vector Vec3; Vec3.x=Vec1.y*Vec2.z-Vec1.z*Vec2.y; Vec3.y=Vec1.z*Vec2.x-Vec1.x*Vec2.z; Vec3.z=Vec1.x*Vec2.y-Vec1.y*Vec2.x; return Vec3;} /************************** crossv ************************/ /************************ dotv ******************/ /* returns the dot product of two vectors, Vec1 and Vec2 */ komplex dotv(vector Vec1, vector Vec2){ komplex dot; dot=(Vec1.x*Vec2.x+Vec1.y*Vec2.y+Vec1.z*Vec2.z); return dot;} /************************** dotv ************************/ /************************ addv ******************/ vector addv(vector Vec1, vector Vec2){ /* Vec1+Vec2 */ vector Vec3; Vec3.x=Vec1.x+Vec2.x; Vec3.y=Vec1.y+Vec2.y; Vec3.z=Vec1.z+Vec2.z; return Vec3;} /************************** addv ************************/ /************************ subv ******************/ vector subv(vector Vec1, vector Vec2){ /* Vec1-Vec2 */ vector Vec3; Vec3.x=Vec1.x-Vec2.x; Vec3.y=Vec1.y-Vec2.y; Vec3.z=Vec1.z-Vec2.z; return Vec3;} /************************** subv ************************/ /************************ conjv ******************/ vector conjv(vector Vec1){/*returns the conjugate of vector Vec1*/ vector Vec3; Vec3.x=conj(Vec1.x); Vec3.y=conj(Vec1.y); Vec3.z=conj(Vec1.z); return Vec3;} /************************** conjv ************************/ /****************************** braket ************************/ /* bra * matrix * ket */ komplex braket(twoby1 bra, twoby1 ket, matrix M1){ twoby1 result1; komplex result2; double pivotket, pivotbra, pivot; pivot = cabs(ket.m); pivotket = cabs(ket.n); if (pivot > pivotket) pivotket = pivot; if (pivotket < LIM_double) pivotket = 1; pivot = cabs(bra.m); pivot = cabs(bra.m); pivotbra = cabs(bra.n); if (pivot > pivotbra) pivotket = pivot; if (pivotbra < LIM_double) pivotbra = 1; ket.m=ket.m/pivotket; ket.n=ket.n/pivotket; bra.m=bra.m/pivotbra; bra.n=bra.n/pivotbra; if (log10(pivotket)+log10(pivotbra) > 290) {pivot = 1.; printf ("overflow in \"braket\" (too thick layers?) #%d\n", pivotoverflow++);} else pivot = pivotket * pivotbra; result1.m = (M1.a*ket.m)+(M1.b*ket.n); result1.n = (M1.c*ket.m)+(M1.d*ket.n); result2 = (bra.m*result1.m)+(bra.n*result1.n); result2 = result2*pivot; return(result2);} /***************************** braket *************************/ /**************************** mtimesv *************************/ /* matrix, M1 times two-by-one vector, ket */ twoby1 mtimesv(matrix M1, twoby1 ket){ twoby1 vec; vec.m = (M1.a*ket.m) + (M1.b*ket.n); vec.n = (M1.c*ket.m) + (M1.d*ket.n); return(vec);} /************************ end mtimesv *************************/ /************************ add2by1 *****************************/ /* adds together two 2by1 vectors */ twoby1 add2by1(twoby1 t1, twoby1 t2){ twoby1 t3; t3.m = t1.m + t2.m; t3.n = t1.n + t2.n; return t3;} /********************** end add2by1 ***************************/ /****************** matrixTIMESsvector ***********************/ svector matrixTIMESsvector (matrix M, svector S) { svector R; R.Vm.m = M.a * S.Vm.m + M.b * S.Vm.n; R.Vm.n = M.c * S.Vm.m + M.d * S.Vm.n; R.Vn.m = M.a * S.Vn.m + M.b * S.Vn.n; R.Vn.n = M.c * S.Vn.m + M.d * S.Vn.n; return R;} /****************** matrixTIMESsvector ***********************/ /****************** smatrixTIMESsvector *********************** supermatrix times supervector (2by1's) svector smatrixTIMESsvector(smatrix S_temp, svector V1){ svector V2; V2.Vm = add2by1(mtimesv(S_temp.Sa, V1.Vm),mtimesv(S_temp.Sb,V1.Vn)); V2.Vn = add2by1(mtimesv(S_temp.Sc, V1.Vm),mtimesv(S_temp.Sd,V1.Vn)); return(V2);} **************** end smatrixTIMESsvector *********************/ /****************** smatrixTIMESsvector *********************** supermatrix times supervector (2by1's) */ svector smatrixTIMESsvector(smatrix S, svector V) {svector E; double re, im, a1, a2, a3, a4, a5, a6, a7, a8, asum; a1 = creal(S.Sa.a) * creal(V.Vm.m); a2 = creal(S.Sa.b) * creal(V.Vm.n); a3 = creal(S.Sb.a) * creal(V.Vn.m); a4 = creal(S.Sb.b) * creal(V.Vn.n); a5 = - cimag(S.Sa.a) * cimag(V.Vm.m); a6 = - cimag(S.Sa.b) * cimag(V.Vm.n); a7 = - cimag(S.Sb.a) * cimag(V.Vn.m); a8 = - cimag(S.Sb.b) * cimag(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); re = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(re) < asum * LIM_double) re = 0.; a1 = creal(S.Sa.a) * cimag(V.Vm.m); a2 = creal(S.Sa.b) * cimag(V.Vm.n); a3 = creal(S.Sb.a) * cimag(V.Vn.m); a4 = creal(S.Sb.b) * cimag(V.Vn.n); a5 = cimag(S.Sa.a) * creal(V.Vm.m); a6 = cimag(S.Sa.b) * creal(V.Vm.n); a7 = cimag(S.Sb.a) * creal(V.Vn.m); a8 = cimag(S.Sb.b) * creal(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); im = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(im) < asum * LIM_double) im = 0.; E.Vm.m = re + I * im; a1 = creal(S.Sa.c) * creal(V.Vm.m); a2 = creal(S.Sa.d) * creal(V.Vm.n); a3 = creal(S.Sb.c) * creal(V.Vn.m); a4 = creal(S.Sb.d) * creal(V.Vn.n); a5 = - cimag(S.Sa.c) * cimag(V.Vm.m); a6 = - cimag(S.Sa.d) * cimag(V.Vm.n); a7 = - cimag(S.Sb.a) * cimag(V.Vn.m); a8 = - cimag(S.Sb.d) * cimag(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); re = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(re) < asum * LIM_double) re = 0.; a1 = creal(S.Sa.c) * cimag(V.Vm.m); a2 = creal(S.Sa.d) * cimag(V.Vm.n); a3 = creal(S.Sb.c) * cimag(V.Vn.m); a4 = creal(S.Sb.d) * cimag(V.Vn.n); a5 = cimag(S.Sa.c) * creal(V.Vm.m); a6 = cimag(S.Sa.d) * creal(V.Vm.n); a7 = cimag(S.Sb.c) * creal(V.Vn.m); a8 = cimag(S.Sb.d) * creal(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); im = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(im) < asum * LIM_double) im = 0.; E.Vm.n = re + I * im; a1 = creal(S.Sc.a) * creal(V.Vm.m); a2 = creal(S.Sc.b) * creal(V.Vm.n); a3 = creal(S.Sd.a) * creal(V.Vn.m); a4 = creal(S.Sd.b) * creal(V.Vn.n); a5 = - cimag(S.Sc.a) * cimag(V.Vm.m); a6 = - cimag(S.Sc.b) * cimag(V.Vm.n); a7 = - cimag(S.Sd.a) * cimag(V.Vn.m); a8 = - cimag(S.Sd.b) * cimag(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); re = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(re) < asum * LIM_double) re = 0.; a1 = creal(S.Sc.a) * cimag(V.Vm.m); a2 = creal(S.Sc.b) * cimag(V.Vm.n); a3 = creal(S.Sd.a) * cimag(V.Vn.m); a4 = creal(S.Sd.b) * cimag(V.Vn.n); a5 = cimag(S.Sc.a) * creal(V.Vm.m); a6 = cimag(S.Sc.b) * creal(V.Vm.n); a7 = cimag(S.Sd.a) * creal(V.Vn.m); a8 = cimag(S.Sd.b) * creal(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); im = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(im) < asum * LIM_double) im = 0.; E.Vn.m = re + I * im; a1 = creal(S.Sc.c) * creal(V.Vm.m); a2 = creal(S.Sc.d) * creal(V.Vm.n); a3 = creal(S.Sd.c) * creal(V.Vn.m); a4 = creal(S.Sd.d) * creal(V.Vn.n); a5 = - cimag(S.Sc.c) * cimag(V.Vm.m); a6 = - cimag(S.Sc.d) * cimag(V.Vm.n); a7 = - cimag(S.Sd.a) * cimag(V.Vn.m); a8 = - cimag(S.Sd.d) * cimag(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); re = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(re) < asum * LIM_double) re = 0.; a1 = creal(S.Sc.c) * cimag(V.Vm.m); a2 = creal(S.Sc.d) * cimag(V.Vm.n); a3 = creal(S.Sd.c) * cimag(V.Vn.m); a4 = creal(S.Sd.d) * cimag(V.Vn.n); a5 = cimag(S.Sc.c) * creal(V.Vm.m); a6 = cimag(S.Sc.d) * creal(V.Vm.n); a7 = cimag(S.Sd.c) * creal(V.Vn.m); a8 = cimag(S.Sd.d) * creal(V.Vn.n); asum=fabs(a1)+fabs(a2)+fabs(a3)+fabs(a4)+fabs(a5)+fabs(a6)+fabs(a7)+fabs(a8); im = a1 + a2 + a3 + a4 + a5 + a6 + a7 + a8; if (fabs(im) < asum * LIM_double) im = 0.; E.Vn.n = re + I * im; return(E);} /**************** end smatrixTIMESsvector *********************/ /****************** smatrixTIMESmatrix ***********************/ smatrix smatrixTIMESmatrix (smatrix S, matrix M) { /*supermatrix times matrix*/ smatrix R; R.Sa = matrixTIMESmatrix(S.Sa, M); R.Sb = matrixTIMESmatrix(S.Sb, M); R.Sc = matrixTIMESmatrix(S.Sc, M); R.Sd = matrixTIMESmatrix(S.Sd, M); return R;} /****************** smatrixTIMESmatrix ***********************/ /****************** matrixTIMESsmatrix ***********************/ smatrix matrixTIMESsmatrix (matrix M, smatrix S) { /*matrix times supermatrix*/ smatrix R; R.Sa = matrixTIMESmatrix(M, S.Sa); R.Sb = matrixTIMESmatrix(M, S.Sb); R.Sc = matrixTIMESmatrix(M, S.Sc); R.Sd = matrixTIMESmatrix(M, S.Sd); return R;} /****************** smatrixTIMESmatrix ***********************/ /****************** smatrixTIMEScomplex ***********************/ smatrix smatrixTIMEScomplex (smatrix S, komplex k){ smatrix R; R.Sa = matrixTIMEScomplex(S.Sa, k); R.Sb = matrixTIMEScomplex(S.Sb, k); R.Sc = matrixTIMEScomplex(S.Sc, k); R.Sd = matrixTIMEScomplex(S.Sd, k); return R;} /****************** smatrixTIMEScomplex ***********************/ /***************************************************************/ /* */ /* Read and Write routines */ /* */ /***************************************************************/ /************************ readminfile ******************/ void readminfile(void) { int nd = 0, ng = 0, counter = 0; char line[MAX_LEN+1], blob[MAX_LEN+1]; float fstart, fstep, fmin, fmax; FILE *file; file=fopen(minfile,"r"); if (file==NULL) terminate(" The minuit command file does not exist! "); while (fgets (line, MAX_LEN, file) != NULL) { if (line[0] == '&') break; if ((line[0] == '#') && (line[1] == '$')) {counter++; if (counter > 1) break;} if (line[0] != '#') { if (line[0] == '$') {ng = ng + 1; sscanf (line, "%s%d%d%s ", blob, &nrep[ng], &pos[ng], layername_in[ng]); if (mode == 9) printf ("%4d %4d %4d \n", ng, nrep[ng], pos[ng]);} else {nd = nd + 1; sscanf(line, "%s%g%g%g%g", parname[nd], &fstart, &fstep, &fmin, &fmax); start[nd] = fstart; step[nd] = fstep; min[nd] = fmin; max[nd] = fmax; if (mode == 9) printf(" %4i %10s %12.4g %12.4g %12.4g %12.4g\n", nd, parname[nd], start[nd], step[nd], min[nd], max[nd]);}}} fclose(file); if (ng*LAY_PARAM+GEN_PARAM != nd) terminate("Format error in the minuit command file !"); else ndlayer = ng; im = nd; /* the number of fit parameters */ printf("The read in of the minuit command file \"%s\" is finished.\n", minfile);} /************************ readminfile ******************/ /************************ readdatfile ******************/ void readdatfile(void) { int i, file_open = 1; char line[MAX_LEN+1]; float f1, f2, f3, f4, f5; FILE *file; file = fopen (datfile, "r"); if(file==NULL) { file_open = 0; if (fit_dim == 0) message ("No data file has been found.",""); else terminate ( "No data file has been found: the fitting is not possible."); if ((mode==2)||(mode==12)) terminate ( "Please use mode #1 or #11 instead of #2 or #12."); } else while (fgets (line,MAX_LEN,file)!=NULL) { if (line[0]!='#') {iexp++; sscanf (line, "%f%f%f%f%f", &f1, &f2, &f3, &f4, &f5); if (f1 < 0) f1 = -f1; alpha_in[iexp] = f1; N[iexp] = f2 * NFACTOR; BG[iexp] = f3 * bg_factor; Err[iexp] = sqrt(fabs(N[iexp]-BG[iexp])); monitor[iexp] = f4; agroup[iexp] = f5; if ((mode == 1)||(mode == 11)) ; else{ if (f4 == 0) terminate ("readdatfile: monitor=0 !"); if (f5 > 5) terminate ("readdatfile: agroup > 5 !"); if (f5 < 1) terminate ("readdatfile: agroup < 1 !");}}} {/* reorganize the dat file in order to obtain a monotonous set of ai: */ int change; double b_alpha, b_N, b_BG, b_monitor, b_Err; int b_agroup; do {change = 0; for (i=1; ialpha_in[i+1]) {change = 1; b_alpha = alpha_in[i]; alpha_in[i] = alpha_in[i+1]; alpha_in[i+1] = b_alpha; b_N = N[i]; N[i]= N[i+1]; N[i+1]=b_N; b_BG = BG[i]; BG[i]= BG[i+1]; BG[i+1]=b_BG; b_Err = Err[i]; Err[i]= Err[i+1]; Err[i+1]=b_Err; b_monitor = monitor[i]; monitor[i]= monitor[i+1]; monitor[i+1]=b_monitor; b_agroup = agroup[i]; agroup[i]= agroup[i+1]; agroup[i+1]=b_agroup;} }} while (change != 0); } if (mode == 9) for (i=1; i<=iexp; i++) printf (" %4i %14.4f %14.4f %14.4f %14.4f %4i\n", i, alpha_in[i], N[i], Err[i], monitor[i], agroup[i]); if (file_open == 1){ printf ("The read in of the data file \"%s\" is finished. \n", datfile); if ((mode == 1)||(mode==11)||(mode==5)||(mode==6)) message ("The data file will not be used in this mode.\n",""); else message("",""); fclose(file);} if (iexp > NE) terminate ("readdatfile: The size of arrays are too small" " for this amount of data: please make NE larger.");} /************************ readdatfile ******************/ /************************ readpoldatfile ******************/ void readpoldatfile(void) { /* This function reads in the results of a four-segment specular polarized neutron reflectivity (PNR) experiment */ int i, file_open = 1; char line[MAX_LEN+1]; float f1, f2, f3, f4, f5, f6, f7, f8, f9, f10; FILE *file; file = fopen (datfile, "r"); if(file==NULL) { file_open = 0; if (fit_dim == 0) message ("No data file has been found.",""); else terminate ( "No data file has been found: the fitting is not possible."); if ((mode==2)||(mode==12)) terminate ( "Please use mode #1 or #11 instead of #2 or #12."); } else while (fgets (line,MAX_LEN,file)!=NULL) { if (line[0]!='#') {iexp++; sscanf (line, "%f%f%f%f%f%f%f%f%f%f", &f1, &f2, &f3, &f4, &f5, &f6, &f7, &f8, &f9, &f10); if (f1 < 0) f1 = -f1; alpha_in[iexp] = f1; Ndu[iexp] = f2 * NFACTOR; Errdu[iexp] = sqrt(fabs(Ndu[iexp])); monitordu[iexp] = f3; Nuu[iexp] = f4 * NFACTOR; Erruu[iexp] = sqrt(fabs(Nuu[iexp])); monitoruu[iexp] = f5; Ndd[iexp] = f6 * NFACTOR; Errdd[iexp] = sqrt(fabs(Ndd[iexp])); monitordd[iexp] = f7; Nud[iexp] = f8 * NFACTOR; Errud[iexp] = sqrt(fabs(Nud[iexp])); monitorud[iexp] = f9; agroup[iexp] = f10; if (f3 == 0) terminate ("readpoldatfile: monitor=0 !"); if (f5 == 0) terminate ("readpoldatfile: monitor=0 !"); if (f7 == 0) terminate ("readpoldatfile: monitor=0 !"); if (f9 == 0) terminate ("readpoldatfile: monitor=0 !"); if (f10 > 5) terminate ("readpoldatfile: agroup > 5 !"); if (f10 < 1) terminate ("readpoldatfile: agroup < 1 !");}} {/* reorganize the dat file in order to obtain a monotonous set of ai: */ int change; double b_alpha, b_Nuu, b_Nud, b_Ndu, b_Ndd, b_monuu, b_monud, b_mondu, b_mondd, b_Erruu, b_Errud, b_Errdu, b_Errdd; int b_agroup; do {change = 0; for (i=1; ialpha_in[i+1]) {change = 1; b_alpha = alpha_in[i]; alpha_in[i] = alpha_in[i+1]; alpha_in[i+1] = b_alpha; b_Nuu = Nuu[i]; Nuu[i]= Nuu[i+1]; Nuu[i+1]=b_Nuu; b_Nud = Nud[i]; Nud[i]= Nud[i+1]; Nud[i+1]=b_Nud; b_Ndu = Ndu[i]; Ndu[i]= Ndu[i+1]; Ndu[i+1]=b_Ndu; b_Ndd = Ndd[i]; Ndd[i]= Ndd[i+1]; Ndd[i+1]=b_Ndd; b_Erruu = Erruu[i]; Erruu[i]= Erruu[i+1]; Erruu[i+1]=b_Erruu; b_Errud = Errud[i]; Errud[i]= Errud[i+1]; Errud[i+1]=b_Errud; b_Errdu = Errdu[i]; Errdu[i]= Errdu[i+1]; Errdu[i+1]=b_Errdu; b_Errdd = Errdd[i]; Errdd[i]= Errdd[i+1]; Errdd[i+1]=b_Errdd; b_monuu = monitoruu[i]; monitoruu[i]= monitoruu[i+1]; monitoruu[i+1]=b_monuu; b_monud = monitorud[i]; monitorud[i]= monitorud[i+1]; monitorud[i+1]=b_monud; b_mondu = monitordu[i]; monitordu[i]= monitordu[i+1]; monitordu[i+1]=b_mondu; b_mondd = monitordd[i]; monitordd[i]= monitordd[i+1]; monitordd[i+1]=b_mondd; b_agroup = agroup[i]; agroup[i]= agroup[i+1]; agroup[i+1]=b_agroup;} }} while (change != 0); } if (file_open == 1){ printf ("The read in of the data file \"%s\" is finished. \n", datfile); fclose(file);} if (iexp > NE) terminate ("readpoldatfile: The size of arrays are too small" " for this amount of data: please make NE larger.");} /************************ readpoldatfile ******************/ /************************* rewrite **************************/ void rewrite (void){ /* after fitting the minuit command file will be rewritten */ char line[MAX_LEN+1], ch[MAX_LEN+1]; int i,j,k; FILE *file; FILE *filb; message(minfilb, "(the copy of the minuit command file) will be written or rewritten."); message(minfile, "will be rewritten."); /* the back-up minuit command file will be written: */ file = fopen(minfile, "r"); filb = fopen(minfilb, "w"); next: if (fgets(line,MAX_LEN,file)!=NULL) {fprintf(filb, "%s", line); goto next;} fclose (filb); fclose (file); /* the minuit command file will be rewritten: */ file = fopen(minfile, "r+"); do {if (fgets(line,MAX_LEN,file)==NULL) terminate("the minfile is too short /n");} while (line[1]!='$'); strcpy(ch, line); fflush(file); for (j=1; j<=GEN_PARAM; j++) { fprintf (file, " %10s %10.4g %10.4g %10.4g %10.4g x[%2i]\n", parname[j], start[j], step[j], min[j], max[j], j-1);} printf("ndlayer: %d\n",ndlayer); for (i=1; i<=ndlayer; i++){ fprintf (file, "$ %d %d %s\n", nrep[i], pos[i], layername_in[i]); for (j=1; j<=LAY_PARAM; j++){ k = GEN_PARAM+((i-1)*LAY_PARAM)+j; fprintf (file, " %10s %10.4g %10.4g %10.4g %10.4g x[%2i]\n", parname[k], start[k], step[k], min[k], max[k], k-1);}} strcat (ch, "#this was the top (reflexion) layer next to vacuum \n"); fprintf(file, "%s", ch); fprintf(file, "#\n"); fprintf(file, "# fcn = %18.2f; fcnred = %10.2f \n", fcn_end, fcnred_end); fprintf(file, "# imin = %5i; imax = %5i; npar = %8i\n", imin, imax, n_par); fprintf(file, "# \n"); for (i=1; i<=10; i++) fprintf(file, " \n"); fclose (file); /* additionally, it will be written to the screen: */ for (j=1; j<=im; j++) printf (" %10s %10.4g %10.4g %10.4g %10.4g x[%2i]\n", parname[j], start[j], step[j], min[j], max[j], j-1); printf(" fcn = %18.2f; fcnred = %10.2f \n", fcn_end, fcnred_end); printf(" imin = %5i; imax = %5i; npar = %8i\n", imin, imax, n_par); message(minfilb, "(the copy of the minuit command file) has been written or rewritten."); message(minfile, "has been rewritten.");} /*************************** rewrite ***************************/ /***************************** simulate ****************************/ void simulate (void) { /* x... (in), y... (final) in degrees */ int jx, jy, kx, ky; double x, y; for (jy = pixel; jy >= 0; jy--) { if(jy%SCREEN_WRITE==0) {printf("%d ",jy); fflush(stdout);} y = y_min + (y_max - y_min) * jy / (double) pixel; for (jx = 0; jx <= pixel; jx++) { x = x_min + (x_max - x_min) * jx / (double) pixel; mod[jx][jy] = model_diffuse(x, y);}} printf("\n"); blur (mod); for (kx=0; kx <= pixel; kx++) for (ky=0; ky <= pixel; ky++) sim[si][kx][ky] = mod[kx][ky];} /************************** simulate ***************************/ /***************************** asymmetry ****************************/ void asymmetry (void) { int jx, jy; for (jy = pixel; jy >= 0; jy--) for (jx = 0; jx <= pixel; jx++) { sim[5][jx][jy] = (sim[1][jx][jy]-sim[0][jx][jy])/(1.e-10+sim[1][jx][jy]+sim[0][jx][jy]); sim[6][jx][jy] = (sim[2][jx][jy]-sim[3][jx][jy])/(1.e-10+sim[2][jx][jy]+sim[3][jx][jy]);} blur (sim[5]); blur (sim[6]); } /***************************** asymmetry ****************************/ /***************************** log_2D ****************************/ void log_2D (void) { int kx, ky; /* The base 10 logarithms of the 2D simulated maps will be calculated: linear scale: sim[][][] -> log scale: lsim[][][] */ /* original version: for (si=0; si<=3; si++) if (si!=4) { for (kx=0; kx <= pixel; kx++) for (ky=0; ky <= pixel; ky++) if (sim[si][kx][ky] > EFFZERO) sim[si][kx][ky] = log10(1.e-5 + sim[si][kx][ky]); else sim[si][kx][ky] = log10(1.e-5);} */ /* test version 15. 07. 2004: */ for (si=0; si<=3; si++) for (kx=0; kx <= pixel; kx++) for (ky=0; ky <= pixel; ky++) lsim[si][kx][ky] = log10(bg_const + sim[si][kx][ky]); } /************************** log_2D ***************************/ /************************** blur ***************************/ /* global parameters: channels: 1 ... pixel in both dimensions blur_i in degree blur_f in degree lower in degree higher in degree */ void blur (double in [MAXPIX+1] [MAXPIX+1]) { double coeff, is_2, js_2; int is, js; int i, j, ii, jj, imin, imax, jmin, jmax; for (i=0; i<=pixel; i++) for (j=0; j<=pixel; j++) { out [i][j] = 0.0; nor [i][j] = 0.0;} is = ceil (pixel*blur_i/(higher-lower)); js = ceil (pixel*blur_f/(higher-lower)); is_2 = 1./(2. * is * is); js_2 = 1./(2. * js * js); for (i=0; i<=pixel; i++) { imin = i - 5*is; if (imin < 1) imin = 1; imax = i + 5*is; if (imax > pixel) imax = pixel; for (j=0; j<=pixel; j++){ jmin = j - 5*js; if (jmin < 1) jmin = 1; jmax = j + 5*js; if (jmax > pixel) jmax = pixel; for (ii=imin; ii<=imax; ii++) for (jj=jmin; jj<=jmax; jj++) { coeff = - (ii-i)*(ii-i)*is_2 - (jj-j)*(jj-j)*js_2; coeff = 4. * log(2.) * coeff; /* that the Gaussian parameter = FWHM is */ coeff = coeff < -LIM_EXP ? 0 : exp(coeff); out [i][j] = out [i][j] + in [ii][jj] * coeff; nor [i][j] = nor [i][j] + coeff;}}} for (i=0; i<=pixel; i++) for (j=0; j<=pixel; j++) if (in [i][j] > 1.e-5) in [i][j] = out [i][j] / nor [i][j]; /* 1e-10 has been checked: no change */ } /************************** blur ***************************/ /************************** blur1D ***************************/ /* This function is probably unnecessary. Presently NOT USED at all. global parameters: channels: 1 ... pixel in one dimension blur_m in degree lower in degree higher in degree */ void blur1D (double in[MAXPIX+1]) { double coeff, is_2; int is; int i, ii, imin, imax; for (i=0; i<=pixel; i++) {out [i][1] = 0.0; nor [i][1] = 0.0;} is = ceil (pixel*blur_m/(higher-lower)); is_2 = 1./(2. * is * is); for (i=0; i<=pixel; i++) { imin = i - 5*is; if (imin < 1) imin = 1; imax = i + 5*is; if (imax > pixel) imax = pixel; for (ii=imin; ii<=imax; ii++) { coeff = - (ii-i)*(ii-i)*is_2; coeff = coeff < -LIM_EXP ? 0 : exp(coeff); out [i][1] = out [i][1] + in [ii] * coeff; nor [i][1] = nor [i][1] + coeff;}} for (i=0; i<=pixel; i++) if (in [i] > 1.e-5) in [i] = out [i][1] / nor [i][1];} /************************** blur1D ***************************/ /****************************** ppm ********************************/ void ppm(void) { FILE *file; double dummy, fdummy; int jx, jy, io, i1; double function_5_min, function_5_max, function_6_min, function_6_max; /* against compiler warning message: */ file = fopen(ppmfile, "w"); fclose( file); si = 0; if ((mode == 5) || (mode == 6)) { /* no fit, only simulation */ for (io=0; io<=1; io++) for (i1=0; i1<=1; i1++) { /*P_polariz = 2 * io - 1; P_analyse = 2 * i1 - 1;*/ if (io == 0) P_polariz = start[2]; if (io == 1) P_polariz = - start[2]; if (i1 == 0) P_analyse = start[4]; if (i1 == 1) P_analyse = - start[4]; N_polariz = start[3]; N_analyse = start[5]; /* = set for simulation purposes only*/ simulate(); si = si + 1; } asymmetry (); log_2D (); /* calculation of the allover function range: */ function_min = lsim[0][1][1]; function_max = lsim[0][1][1]; for (si=0; si<=3; si++) for (jy = pixel; jy >= 1; jy--) for (jx = 1; jx <= pixel; jx++) { if (function_min > lsim[si][jx][jy]) function_min = lsim[si][jx][jy]; if (function_max < lsim[si][jx][jy]) function_max = lsim[si][jx][jy];} 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"); function_5_min = sim[5][1][1]; function_5_max = sim[5][1][1]; function_6_min = sim[6][1][1]; function_6_max = sim[6][1][1]; for (jy = pixel; jy >= 1; jy--) for (jx = 1; jx <= pixel; jx++) { if (function_5_min > sim[5][jx][jy]) function_5_min = sim[5][jx][jy]; if (function_5_max < sim[5][jx][jy]) function_5_max = sim[5][jx][jy]; if (function_6_min > sim[6][jx][jy]) function_6_min = sim[6][jx][jy]; if (function_6_max < sim[6][jx][jy]) function_6_max = sim[6][jx][jy];} fprintf(stderr, "au function range: min = %g; max = %g \n", function_5_min, function_5_max); fprintf(stderr, "ad function range: min = %g; max = %g \n", function_6_min, function_6_max); for (jy = pixel; jy >= 1; jy--) for (jx = 1; jx <= pixel; jx++) sim[4][jx][jy] = function_min + /* (pixel - jx) * (function_max - function_min )/(double) pixel; */ (jx) * (function_max - function_min )/(double) pixel; for (si = 0; si <=6; si++){ if (si == 0) file = fopen(ppmfile_uu, "w"); if (si == 1) file = fopen(ppmfile_ud, "w"); if (si == 2) file = fopen(ppmfile_du, "w"); if (si == 3) file = fopen(ppmfile_dd, "w"); if (si == 4) file = fopen(ppmfile, "w"); if (si == 5) file = fopen(ppmfile_au, "w"); if (si == 6) file = fopen(ppmfile_ad, "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"); if (si == 5) fdummy = (colour_max - colour_min)/(function_5_max-function_5_min); else if (si == 6) fdummy = (colour_max - colour_min)/(function_6_max-function_6_min); else fdummy = (colour_max - colour_min)/(function_max-function_min); for (jy = pixel; jy >= 1; jy--) for (jx = 1; jx <= pixel; jx++) { /* pixel definitions: */ if (si == 5) dummy = (sim[si][jx][jy] - function_5_min) * fdummy + colour_min; else if (si == 6) /* itt eredetileg 5 allt a 6 helyett:*/ dummy = (sim[si][jx][jy] - function_6_min) * fdummy + colour_min; else if (si == 4) dummy = (sim[si][jx][jy] - function_min) * fdummy + colour_min; else dummy = (lsim[si][jx][jy] - function_min) * fdummy + colour_min; colour (file, dummy, 1., 1.); } /* closing with \n: */ fprintf (file, " \n"); fclose(file);} /* write ascii file: */ for (si = 0; si <=3; si++){ if (si == 0) file = fopen(ascfile_uu, "w"); if (si == 1) file = fopen(ascfile_ud, "w"); if (si == 2) file = fopen(ascfile_du, "w"); if (si == 3) file = fopen(ascfile_dd, "w"); for (jy = pixel; jy >= 1; jy--){ for (jx = 1; jx <= pixel; jx++) { fprintf (file, "%g ", sim[si][jx][jy]);} fprintf (file, "\n");} fclose(file);} } printf ("dwba overflow = %10d times; dwba ok = %10d times\n", dwba_error, dwba_ok); /* fprintf(stderr, "The size of the picture in %s is %i x %i pixel \n", ppmfile, pixel, pixel); */ message("",""); message(ppmfile_uu, "has been written or rewritten,"); message(ppmfile_ud, "has been written or rewritten,"); message(ppmfile_du, "has been written or rewritten,"); message(ppmfile_dd, "has been written or rewritten,"); message(ascfile_uu, "has been written or rewritten,"); message(ascfile_ud, "has been written or rewritten,"); message(ascfile_du, "has been written or rewritten,"); message(ascfile_dd, "has been written or rewritten,"); message(ppmfile_au, "has been written or rewritten,"); message(ppmfile_ad, "has been written or rewritten,"); message(ppmfile, " has been written or rewritten."); message("",""); message("You may want to run the 'gqview' program", "to display these 2D plots.");} /*************************** 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=floor(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)(floor(255.*r)), (int)(floor(255.*g)), (int)(floor(255.*b))); } /**************************** colour *******************************/ /*********************** grafika **********************/ /* GNUPLOT HOWTO: t '' w l (in newer gnuplot versions:) lw lt (in older gnuplot versions:) lw w e == with errorbars */ void grafika(void){ int j, io, i1, point; double e[5][2000]; double buffer; char ch[MAX_LEN], dum; double von, bis, von1, bis1, delta; /* zoom parameters */ int zoom=1; /* zoom factor */ double zn=1; /* zoom range */ FILE *file; FILE *pipe; model_specular_no_blurring (); if ((mode == 1)||(mode == 2)) {/* no fit, only simulation, 4 segments */ for (io=0; io<=1; io++) { for (i1=0; i1<=1; i1++) { /* P_polariz = 2 * io - 1; P_analyse = 2 * i1 - 1; N_polariz = 0.5; N_analyse = 0.5; */ if (io == 0) P_polariz = start[2]; if (io == 1) P_polariz = - start[2]; if (i1 == 0) P_analyse = start[4]; if (i1 == 1) P_analyse = - start[4]; N_polariz = start[3]; N_analyse = start[5]; /* = set for simulation purposes only*/ model_specular_no_blurring (); for (j=0; j<=pixel; j++) { buffer = lower + j*(higher-lower)/(double) pixel; e[1+io+2*i1][j]=model_specular(buffer); e[0][j]=buffer;}}} file = fopen(simfile, "w"); for (j=0; j<=pixel; j++) fprintf (file, "%+-16.8e %+-16.8e %+-16.8e %+-16.8e %+-16.8e\n", e[0][j], e[1][j], e[2][j], e[3][j], e[4][j]); /* -- -+ +- ++ */ fclose(file);} else if ((mode == 11)||(mode == 12)) {/* no fit, only simulation */ /* P_polariz = start [2]; P_analyse = start [4]; N_polariz = start [3]; N_analyse = start [5]; these have been initialized anyhow */ model_specular_no_blurring (); for (j=0; j<=pixel; j++) { buffer = lower + j*(higher-lower)/(double) pixel; e[1][j]=model_specular(buffer); e[0][j]=buffer;} file = fopen(simfile, "w"); for (j=0; j<=pixel; j++) fprintf (file, "%+-16.8e %+-16.8e \n", e[0][j], e[1][j]); fclose(file);} else { /* fit and experimental points */ file = fopen(simfile, "w"); point = pixel * (imax - imin)/iexp; for (j=1; j<=point; j++) { buffer = alpha[imin]+j*(alpha[imax]-alpha[imin])/(double) point; fprintf (file, "%g \t %g \n", buffer, model_specular(buffer));} fclose(file);} file = fopen(modfile, "w"); for (j=0; j<=pixel; j++) { buffer = alpha[1]+j*(alpha[iexp]-alpha[1])/(double) pixel; if ((buffer<=alpha[imin])||(buffer>=alpha[imax])) fprintf (file, "%g \t %g \n", buffer, model_specular(buffer)); else fprintf (file, "\n");} fclose(file); file = fopen(expfile, "w"); for (j=1; j<=iexp; j++) fprintf (file, "%16.4g %16.4g %16.4g \n", alpha[j], (N[j]-BG[j])/(monitor[j]*absorber[agroup[j]]), Err[j]/(monitor[j]*absorber[agroup[j]])); fclose(file); nbplot (); printf ("The file \"%s\" has been written or rewritten. \n", nbfile); do {pipe = popen("gnuplot -geometry 900x400+10+10", "w"); fprintf(pipe, "set grid \n"); fflush(pipe); if (ch[0] == 'p') { fprintf(pipe, "set terminal postscript color solid lw 2\n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", psfile); fflush(pipe);} if (ch[0] == 'c') { fprintf(pipe, "set terminal postscript color solid lw 2\n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", nbpsfile); fflush(pipe);} if (ch[0] == 'g') {if (gnuplotversion == 0) {fprintf(pipe, "set terminal png \n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", pngfile); fflush(pipe);} else {};} if (ch[0] == 'b') {if (gnuplotversion == 0) {fprintf(pipe, "set terminal png \n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", nbpngfile); fflush(pipe);} else {};} if ((mode==1)||(mode==2)||(mode==11)||(mode==12)) fprintf (pipe, "set title 'Reflectivity <%s> [counts/monitor counts vs. angle in degrees]' \n", name); else fprintf (pipe, "set title 'minuit fit [%s]' \n", name); fflush(pipe); fprintf (pipe, "set logscale y \n"); fflush(pipe); if ((ch[0]== 'n') || (ch[0]== 'c') || (ch[0] == 'b')){ /* scattering length density plot: */ fprintf (pipe, "set title 'Scattering-length density vs. depth <%s> [A^-2 vs. A]' \n", name); fprintf (pipe, "set nologscale y \n"); fflush(pipe); fprintf (pipe, "set xlabel 'reflection surface " " substrate >' \n"); fflush(pipe); if (gnuplotversion == 0) fprintf (pipe, "plot [:][:] '%s' t '' w l lw 2 lt 1 \n", nbfile); else fprintf (pipe, "plot [:][:] '%s' t '' w l lw 2 1 \n", nbfile); fflush(pipe);} else {/* reflectivity plot: */ /* zoom calculation: */ if (ch[0] == 'z') {sscanf (ch, "%1c%1i", &dum, &zoom); if (zoom < 1) zoom = 1; zn = (zoom+1.)/2.;} if (ch[0] == '+') zn = zn + .5; else if (ch[0] == '=') zn = zn + .5; else if (ch[0] == '-') zn = zn - .5; if (zn < 1) zn = 1; if (zn > zoom) zn = zoom; von1 = lower; bis1 = higher; delta = (bis1 - von1)/zoom; von = von1 + (zn-1) * delta; bis = von + delta; if (mode==1) {if (gnuplotversion == 0) fprintf (pipe, "plot [%f:%f][:] '%s' using ($1):($2) t '--' w l lw 2 lt 1, " "'%s' using ($1):($3) t '-+' w l lw 2 lt 3, " "'%s' using ($1):($4) t '+-' w l lw 2 lt 2, " "'%s' using ($1):($5) t '++' w l lw 2 lt 7 \n", von, bis, simfile, simfile, simfile, simfile); else fprintf (pipe, "plot [%f:%f][:] '%s' using ($1):($2) t '--' w l lw 2 1, " "'%s' using ($1):($3) t '-+' w l lw 2 3, " "'%s' using ($1):($4) t '+-' w l lw 2 2, " "'%s' using ($1):($5) t '++' w l lw 2 7 \n", von, bis, simfile, simfile, simfile, simfile);} else if (mode==11) {if (gnuplotversion == 0) fprintf (pipe, "plot [%f:%f][:] '%s' using ($1):($2) t '' w l lw 2 lt 1, " "'%s' using ($1):($3) t '' w l lw 2 lt 3, " "'%s' using ($1):($4) t '' w l lw 2 lt 2, " "'%s' using ($1):($5) t '' w l lw 2 lt 7 \n", von, bis, simfile, simfile, simfile, simfile); else fprintf (pipe, "plot [%f:%f][:] '%s' using ($1):($2) t '' w l lw 2 1, " "'%s' using ($1):($3) t '' w l lw 2 3, " "'%s' using ($1):($4) t '' w l lw 2 2, " "'%s' using ($1):($5) t '' w l lw 2 7 \n", von, bis, simfile, simfile, simfile, simfile);} else if (mode==2) {if (gnuplotversion == 0) fprintf (pipe, "plot [%f:%f][:] '%s' notitle w e lt 3, " "'%s' using ($1):($2) t '--' w l lw 2 lt 1, " "'%s' using ($1):($3) t '-+' w l lw 2 lt 3, " "'%s' using ($1):($4) t '+-' w l lw 2 lt 2, " "'%s' using ($1):($5) t '++' w l lw 2 lt 7 \n", von, bis, expfile, simfile, simfile, simfile, simfile); else fprintf (pipe, "plot [%f:%f][:] '%s' notitle w e 3, " "'%s' using ($1):($2) t '--' w l lw 2 1, " "'%s' using ($1):($3) t '-+' w l lw 2 3, " "'%s' using ($1):($4) t '+-' w l lw 2 2, " "'%s' using ($1):($5) t '++' w l lw 2 7 \n", von, bis, expfile, simfile, simfile, simfile, simfile);} else if (mode==12) {if (gnuplotversion == 0) fprintf (pipe, "plot [%f:%f][:] '%s' notitle w e lt 3, " "'%s' using ($1):($2) t '' w l lw 2 lt 1, " "'%s' using ($1):($3) t '' w l lw 2 lt 3, " "'%s' using ($1):($4) t '' w l lw 2 lt 2, " "'%s' using ($1):($5) t '' w l lw 2 lt 7 \n", von, bis, expfile, simfile, simfile, simfile, simfile); else fprintf (pipe, "plot [%f:%f][:] '%s' notitle w e 3, " "'%s' using ($1):($2) t '' w l lw 2 1, " "'%s' using ($1):($3) t '' w l lw 2 3, " "'%s' using ($1):($4) t '' w l lw 2 2, " "'%s' using ($1):($5) t '' w l lw 2 7 \n", von, bis, expfile, simfile, simfile, simfile, simfile);} else {if (gnuplotversion == 0) fprintf (pipe, "plot [%f:%f][:] '%s' notitle w e lt 3, " "'%s' notitle w l lw 2 lt 1, " "'%s' notitle w l lw 2 lt 2 \n", von, bis, expfile, simfile, modfile); else fprintf (pipe, "plot [%f:%f][:] '%s' notitle w e 3, " "'%s' notitle w l lw 2 1, " "'%s' notitle w l lw 2 2 \n", von, bis, expfile, simfile, modfile);} fflush(pipe);} if (ch[0] == 'p') { printf ("The file \"%s\" has been written or rewritten. \n", psfile); ch[0] = 'b'; goto ende;} if (ch[0] == 'b') { if (gnuplotversion == 0) printf ("The file \"%s\" has been written or rewritten. \n", nbpngfile); ch[0] = 'g'; goto ende;} if (ch[0] == 'g') { if (gnuplotversion == 0) printf ("The file \"%s\" has been written or rewritten. \n", pngfile); ch[0] = 'c'; goto ende;} if (ch[0] == 'c') { printf ("The file \"%s\" has been written or rewritten. \n", nbpsfile); ch[0] = 'q'; goto ende;} if ((ch[0] != 'p')|| (ch[0] != 'g') || (ch[0] != 'b') || (ch[0] != 'c')) {printf("'RETURN' (= reflectivity plot) \ 'p' (= makes ps and png files, quit)\n" "'n' (= nb plot) 'q' (= quit) 'z4'|'+'|'-' (= zoom) --->"); fgets (ch, MAX_LEN, stdin);} ende: fprintf (pipe, "q \n"); fflush(pipe); } while (ch[0] != 'q'); fclose (pipe);} /************************* grafika ***************************/ /*********************** grafikaPNR **********************/ void grafikaPNR(void){ int j, io, i1; double e[5][2000]; double buffer, P_pol0, P_ana0; char ch[MAX_LEN], dum; double von, bis, von1, bis1, delta; /* zoom parameters */ int zoom=1; /* zoom factor */ double zn=1; /* zoom range */ FILE *file; FILE *pipe; model_specular_no_blurring (); /* PNR fit, 4 segments */ P_pol0 = P_polariz; P_ana0 = P_analyse; for (io=0; io<=1; io++) { for (i1=0; i1<=1; i1++) { if (io == 0) P_polariz = P_pol0; if (io == 1) P_polariz =-P_pol0; if (i1 == 0) P_analyse = P_ana0; if (i1 == 1) P_analyse =-P_ana0; model_specular_no_blurring (); for (j=0; j<=pixel; j++) { buffer = lower + j*(higher-lower)/(double) pixel; e[1+io+2*i1][j]=model_specular(buffer); e[0][j]=buffer;} } } file = fopen(simfile, "w"); for (j=0; j<=pixel; j++) fprintf (file, "%+-16.8e %+-16.8e %+-16.8e %+-16.8e %+-16.8e\n", e[0][j], e[2][j], e[1][j], e[4][j], e[3][j]); /* du uu dd ud */ fclose(file); P_polariz = P_pol0 ; P_analyse = P_ana0; /* ide jon az adat plot file: */ file = fopen(expfile, "w"); for (j=1; j<=iexp; j++) fprintf (file,"%16.4g %16.4g %16.4g %16.4g %16.4g" "%16.4g %16.4g %16.4g %16.4g \n", alpha[j], Ndu[j]/(monitordu[j]*absorber[agroup[j]]), Errdu[j]/(monitordu[j]*absorber[agroup[j]]), Nuu[j]/(monitoruu[j]*absorber[agroup[j]]), Erruu[j]/(monitoruu[j]*absorber[agroup[j]]), Ndd[j]/(monitordd[j]*absorber[agroup[j]]), Errdd[j]/(monitordd[j]*absorber[agroup[j]]), Nud[j]/(monitorud[j]*absorber[agroup[j]]), Errud[j]/(monitorud[j]*absorber[agroup[j]])); fclose(file); // file = fopen(modfile, "w"); // for (j=0; j<=pixel; j++) { // buffer = alpha[1]+j*(alpha[iexp]-alpha[1])/(double) pixel; // if ((buffer<=alpha[imin])||(buffer>=alpha[imax])) // fprintf (file, "%g \t %g \n", buffer, model_specular(buffer)); // else fprintf (file, "\n");} // fclose(file); nbplot (); printf ("The file \"%s\" has been written or rewritten. \n", nbfile); do {pipe = popen("gnuplot -geometry 900x400+10+10", "w"); fprintf(pipe, "set grid \n"); fflush(pipe); if (ch[0] == 'p') { fprintf(pipe, "set terminal postscript color solid lw 2\n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", psfile); fflush(pipe);} if (ch[0] == 'c') { fprintf(pipe, "set terminal postscript color solid lw 2\n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", nbpsfile); fflush(pipe);} if (ch[0] == 'g') {if (gnuplotversion == 0) {fprintf(pipe, "set terminal png \n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", pngfile); fflush(pipe);} else {};} if (ch[0] == 'b') {if (gnuplotversion == 0) {fprintf(pipe, "set terminal png \n"); fflush(pipe); fprintf(pipe, "set output '%s' \n", nbpngfile); fflush(pipe);} else {};} fprintf (pipe, "set title 'minuit fit [%s]' \n", name); fflush(pipe); fprintf (pipe, "set logscale y \n"); fflush(pipe); if ((ch[0]== 'n') || (ch[0]== 'c') || (ch[0] == 'b')){ /* scattering length density plot: */ fprintf (pipe, "set title 'Scattering-length density vs. depth <%s> [A^-2 vs. A]' \n", name); fprintf (pipe, "set nologscale y \n"); fflush(pipe); fprintf (pipe, "set xlabel 'reflection surface " " substrate >' \n"); fflush(pipe); if (gnuplotversion == 0) fprintf (pipe, "plot [:][:] '%s' t '' w l lw 2 lt 1 \n", nbfile); else fprintf (pipe, "plot [:][:] '%s' t '' w l lw 2 1 \n", nbfile); fflush(pipe);} else {/* reflectivity plot: */ /* zoom calculation: */ if (ch[0] == 'z') {sscanf (ch, "%1c%1i", &dum, &zoom); if (zoom < 1) zoom = 1; zn = (zoom+1.)/2.;} if (ch[0] == '+') zn = zn + .5; else if (ch[0] == '=') zn = zn + .5; else if (ch[0] == '-') zn = zn - .5; if (zn < 1) zn = 1; if (zn > zoom) zn = zoom; von1 = lower; bis1 = higher; delta = (bis1 - von1)/zoom; von = von1 + (zn-1) * delta; bis = von + delta; {if (gnuplotversion == 0) fprintf (pipe, "plot [%f:%f][:]" "'%s' using ($1):($2) t 'du' lw 2 lt 1, " "'%s' using ($1):($4) t 'uu' lw 2 lt 3, " "'%s' using ($1):($6) t 'dd' lw 2 lt 2, " "'%s' using ($1):($8) t 'ud' lw 2 lt 7, " "'%s' using ($1):($2) t 'du' w l lw 2 lt 1, " "'%s' using ($1):($3) t 'uu' w l lw 2 lt 3, " "'%s' using ($1):($4) t 'dd' w l lw 2 lt 2, " "'%s' using ($1):($5) t 'ud' w l lw 2 lt 7 \n", von, bis, expfile, expfile, expfile, expfile, simfile, simfile, simfile, simfile); else fprintf (pipe, "plot [%f:%f][:]" "'%s' using ($1):($2) t 'du' 1, " "'%s' using ($1):($4) t 'uu' 3, " "'%s' using ($1):($6) t 'dd' 2, " "'%s' using ($1):($8) t 'ud' 7, " "'%s' using ($1):($2) t 'du' w l lw 2 1, " "'%s' using ($1):($3) t 'uu' w l lw 2 3, " "'%s' using ($1):($4) t 'dd' w l lw 2 2, " "'%s' using ($1):($5) t 'ud' w l lw 2 7 \n", von, bis, expfile, expfile, expfile, expfile, simfile, simfile, simfile, simfile);} fflush(pipe);} if (ch[0] == 'p') { printf ("The file \"%s\" has been written or rewritten. \n", psfile); ch[0] = 'b'; goto ende;} if (ch[0] == 'b') { if (gnuplotversion == 0) printf ("The file \"%s\" has been written or rewritten. \n", nbpngfile); ch[0] = 'g'; goto ende;} if (ch[0] == 'g') { if (gnuplotversion == 0) printf ("The file \"%s\" has been written or rewritten. \n", pngfile); ch[0] = 'c'; goto ende;} if (ch[0] == 'c') { printf ("The file \"%s\" has been written or rewritten. \n", nbpsfile); ch[0] = 'q'; goto ende;} if ((ch[0] != 'p')|| (ch[0] != 'g') || (ch[0] != 'b') || (ch[0] != 'c')) {printf("'RETURN' (= reflectivity plot) \ 'p' (= makes ps and png files, quit)\n" "'n' (= nb plot) 'q' (= quit) 'z4'|'+'|'-' (= zoom) --->"); fgets (ch, MAX_LEN, stdin);} ende: fprintf (pipe, "q \n"); fflush(pipe); } while (ch[0] != 'q'); fclose (pipe);} /************************* grafikaPNR ***************************/ /*********************** nbplot **********************/ void nbplot (void){ int j, i; double DS = 0., D, max, delta, d0, x, NB, SIGMA, dummy, D0; FILE *file; for (j = 1; j <= alllayer; j++) DS = DS + d[j]; D = DS * (1 + BORDER/100.); D0 = DS * (1 + BORDER/200.); d0 = DS * BORDER/200.; max = d0 + d[1]/2.; delta = D / pixel; SIGMA = rough[0]; i = 0; file = fopen(nbfile, "w"); for (j = 0; j <= pixel; j++){x = j * delta; if (x > max) { if (i == alllayer-1) { SIGMA = rough[i+1]; d0 = d0 + d[i+1]; max = D; i = i + 1;} else { SIGMA = rough[i+1]; d0 = d0 + d[i+1]; max = max + (d[i+1]+d[i+2])/2.; i = i + 1;}} if (SIGMA * LIM_EXP < fabs(x-d0)) if ((x - d0) < 0) dummy = - 1.; else dummy = + 1.; else dummy = erf ((x - d0) / SIGMA); NB = (nb[i+1] + nb[i])/2. + (nb[i+1] - nb[i])/2. * dummy; NB = NB + creal(nb_vacuum); fprintf (file, "%+-16.8e %+-16.8e \n", (D0-x) * 1.e10, NB * 1.e-20);} fflush(file); fclose(file);} /************************* nbplot ***************************/ /************************* to_grad ***************************/ double to_grad(double x){return x*180./PI;} /************************* to_grad ***************************/ /************************* to_radian ***************************/ double to_radian(double x){return x*PI/180.;} /************************* to_radian ***************************/ /************************* alpha_to_qz ***************************/ double alpha_to_qz(double alpha_i) /*in error x2*/ {return (2.*PI/wavelen)*sin(to_radian(alpha_i));} /************************* alpha_to_qz ***************************/ /************************* qz_to_alpha ***************************/ double qz_to_alpha(double qz) /*in error x2*/ {return to_grad(asin(qz*wavelen/2./PI));} /************************* qz_to_alpha ***************************/ /************************* rcsin ***************************/ komplex rcsin(komplex z){ /* calculates sin(a+ib)/cosh(b) */ double a, b; a = creal(z); /* a = a - 2*PI* (floor (a/(2*PI))); */ b = cimag(z); return /* cosh(b) */ (sin(a) + I * cos(a) * tanh(b));} /************************* rcsin ***************************/ /************************* rccos ***************************/ komplex rccos(komplex z){ /* calculates cos(a+ib)/cosh(b) */ double a, b; a = creal(z); /* a = a - 2*PI* (floor (a/(2*PI))); */ b = cimag(z); return /* cosh(b) */ (cos(a) - I * sin(a) * tanh(b));} /************************* rccos ***************************/ /************************* ksqrt ********************* komplex ksqrt (komplex x){ komplex z; double xi, xr, zi, zr; xi = cimag (x); xr = creal (x); if (xi == 0) { if (xr < 0) {zr = 0.; zi = sqrt(-xr);} else {zr = sqrt(xr); zi = 0.;}} else {zr = sqrt(0.5*(xr+sqrt(xr*xr+xi*xi))); zi = xi/(2.*zr);} z = zr + I * zi; return z;} ************************* ksqrt *********************/ /************************ apvector ************************/ void apvector (char *text, vector x){ printf("vector: %s = %10.2e %10.2e %10.2e \n", text, cabs(x.x), cabs(x.y), cabs(x.z));} void apsmatrix (char *text, smatrix x){double e=0.; /* double y1, y2, y3, y4; y1 = cabs(x.Sa.a)+cabs(x.Sa.b)+cabs(x.Sa.c)+cabs(x.Sa.d); y2 = cabs(x.Sb.a)+cabs(x.Sb.b)+cabs(x.Sb.c)+cabs(x.Sb.d); y3 = cabs(x.Sc.a)+cabs(x.Sc.b)+cabs(x.Sc.c)+cabs(x.Sc.d); y4 = cabs(x.Sd.a)+cabs(x.Sd.b)+cabs(x.Sd.c)+cabs(x.Sd.d); printf("smatrix: %s = %10.2e %10.2e %10.2e %10.2e \n", text, y1, y2, y3, y4);*/ printf ("S-matrix: (%s)\n", text); printf ("Sa: %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei \n", creal(x.Sa.a)-e, cimag(x.Sa.a)-e, creal(x.Sa.b)-e, cimag(x.Sa.b)-e, creal(x.Sa.c)-e, cimag(x.Sa.c)-e, creal(x.Sa.d)-e, cimag(x.Sa.d)-e); printf ("Sb: %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei \n", creal(x.Sb.a)-e, cimag(x.Sb.a)-e, creal(x.Sb.b)-e, cimag(x.Sb.b)-e, creal(x.Sb.c)-e, cimag(x.Sb.c)-e, creal(x.Sb.d)-e, cimag(x.Sb.d)-e); printf ("Sc: %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei \n", creal(x.Sc.a)-e, cimag(x.Sc.a)-e, creal(x.Sc.b)-e, cimag(x.Sc.b)-e, creal(x.Sc.c)-e, cimag(x.Sc.c)-e, creal(x.Sc.d)-e, cimag(x.Sc.d)-e); printf ("Sd: %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei %+.2e+%.2ei \n", creal(x.Sd.a)-e, cimag(x.Sd.a)-e, creal(x.Sd.b)-e, cimag(x.Sd.b)-e, creal(x.Sd.c)-e, cimag(x.Sd.c)-e, creal(x.Sd.d)-e, cimag(x.Sd.d)-e);} /************************ apvector ************************/ /************************ apmatrix ************************/ void apmatrix (char *text, matrix x){ double y1, y2, y3, y4; y1 = cabs(x.a); y2 = cabs(x.b); y3 = cabs(x.c); y4 = cabs(x.d); printf("matrix: %s = %10.2e %10.2e %10.2e %10.2e \n", text, y1, y2, y3, y4);} /************************ apmatrix ************************/ /************************ aptwoby1 ************************/ void aptwoby1 (char *text, twoby1 x){ double y1, y2; y1 = cabs(x.m); y2 = cabs(x.n); printf("matrix: %s = %10.2e %10.2e \n", text, y1, y2);} /************************ aptwoby1 ************************/ /************************ test ************************/ void test (void) { /* all four possibilities will be calculated */ double refl; int i, io, i1; double control; message("test mode:\n","\n"); for (io=0; io<=00; io++) for (i1=1; i1<=01; i1++) { P_polariz = 2 * io - 1; P_analyse = 2 * i1 - 1; printf("pol-anal: %.4f %.4f %.4f %.4f \n", P_polariz, N_polariz, P_analyse, N_analyse); refl = defSpec(lower); printf ("spin=%1d%1d in=out=%4.2f R_spek=%8.2e \n", io, i1, lower, refl); refl = defSpec(higher); printf ("spin=%1d%1d in=out=%4.2f R_spek=%8.2e \n", io, i1, higher, refl); refl = defOffSpec(lower, higher); printf ("spin=%1d%1d in=%4.2f out=%4.2f R_diff=%8.2e \n", io, i1, lower, higher, refl); printf ("layer=%3d: vacuum \n", alllayer + 1); for (i=alllayer; i>0; i--) { control = cabs(detm(matrixSUBmatrix(matrixTIMESmatrix(Sm[i].Sa,Sm[i].Sd),matrixTIMESmatrix(Sm[i].Sb,Sm[i].Sc)))); if (i <= alllayer) {/* apsmatrix("Sm ", Sm[i]); */ printf(" sm.a = %.8g det_sm = %.8g det_Sm = %.8g\n", cabs(sm[i].a), cabs(detm(sm[i])), cabs(detm(matrixSUBmatrix(matrixTIMESmatrix(Sm[i].Sa,Sm[i].Sd),matrixTIMESmatrix(Sm[i].Sb,Sm[i].Sc)))) );} printf("(%2d-%2d) inup%+.2e%+.2ei %+.2e%+.2ei [%.4e] \n", i, i-1, creal(chicombi_in[i].Vm.m), cimag(chicombi_in[i].Vm.m), creal(chicombi_in[i].Vn.m), cimag(chicombi_in[i].Vn.m), cabs(chicombi_in[i].Vm.m)); printf("(%2d-%2d) indo%+.2e%+.2ei %+.2e%+.2ei [%.4e] \n", i, i-1, creal(chicombi_in[i].Vm.n), cimag(chicombi_in[i].Vm.n), creal(chicombi_in[i].Vn.n), cimag(chicombi_in[i].Vn.n), cabs(chicombi_in[i].Vm.n)); printf("(%2d-%2d) ouup%+.2e%+.2ei %+.2e%+.2ei [%.4e] \n", i, i-1, creal(chicombi_out[i].Vm.m), cimag(chicombi_out[i].Vm.m), creal(chicombi_out[i].Vn.m), cimag(chicombi_out[i].Vn.m), cabs(chicombi_out[i].Vm.m)); printf("(%2d-%2d) oudo%+.2e%+.2ei %+.2e%+.2ei [%.4e] \n", i, i-1, creal(chicombi_out[i].Vm.n), cimag(chicombi_out[i].Vm.n), creal(chicombi_out[i].Vn.n), cimag(chicombi_out[i].Vn.n), cabs(chicombi_in[i].Vm.n)); } printf ("layer= 0: substrat \n"); apsmatrix("Stot", Stot); printf (" detStot = %.4g\n\n", /*cabs(detm(matrixSUBmatrix(matrixTIMESmatrix(Stot.Sa,Stot.Sd),matrixTIMESmatrix(Stot.Sb,Stot.Sc)))),detStot);*/ detStot );}} /************************ test ************************/ /********************** create sample files ******************/ void create_sample_files (){ FILE *file; file = fopen("sample.min", "w"); fprintf(file, "# \n" "# \n" "# This is the 'sample.min' example minuit command file \n" "# for the program 'superfit'. \n" "# This file describes a non-real layered system of more than 160 layers;\n" "# for a 2D demonstration (off-specular simulation) the \n" "# superfit sample 5 200 0 6 \n" "# can be run; the execution may last for 1 hour or shorter. \n" "# For a fast demonstration the \n" "# superfit sample 5 50 0 6 \n" "# will last 4 minutes or shorter, however the picture will be too small. \n" "# For a 1D demonstration (specular simulation) execute \n" "# superfit sample 1 500 0 6 \n" "# which may last 10 seconds or shorter. \n" "# The execution times have been measured on a 700 MHz laptop. \n" "# \n" "# \n" "# \n" "# \n" "# \n" "# The last column (x[]) in the data table will be added or corrected \n" "# automatically by the 'superfit', if the minuit data file is rewritten; \n" "# it is unnecessary to edit. \n" "# \n" "# \n" "# \n" "# The first 30 lines describe the general parameters and the substrat.\n" "# The next layer contacts the substrat.\n" "# The last layer in the list is the top layer bordering the vacuum. \n" "# \n" "# $ nrep pos: nrep = number of repetitions of the group \n" "# pos = number within the group \n" "# \n" "# The background is calculated in the function untergrund(): \n" "# BACKGROUND = bg _const + bg_ampl * exp(-alpha_i/bg_alpha_i)\n" "# \n" "# Units: lambda: meter (wavelength) \n" "# P_polariz: (polarizer efficiency) \n" "# N_polariz: (polarizer intensity gain) \n" "# P_analize: (analyser efficiency) \n" "# N_analize: (analyser intensity gain) \n" "# blur_i, _f: degree (angular blurring width) \n" "# chsi_x,-y,-z: meter (roughness correlation lengths) \n" "# ai_offset degree (if the experimental ai zero is wrong) \n" "# bg_factor (0/1: no/full experimental background) \n" "# nb: meter^-2 (scattering length density) \n" "# nb_vac: it will only be used in the plot part nbplot (real only) \n" "# for the fit or for the 2D staff nb_vac = 0 \n" "# roughness of the substrat: meter \n" "# roughness: in unit of the layer thickness \n" "# d: meter (layer thickness) \n" "# B: tesla (magnitude of the magnetic field) \n" "# fi: degree (for the guide field (y): fi=0 th=0) \n" "# th: degree (for in-plane fields (x-y): th=0) \n" "# alpha_foot: degree (footprint limit angle) \n" "# bg_ampl: count (decaying contribution to background) \n" "# bg_alpha_i: degree (exponential decay constant of the background\n" "# bg_const: count (constant contribution to background \n" "# w2harmonic: (0 ... 1: weight of the 2nd harmonics in beam)\n" "# absorber1 .. 5 multiplies the corresponding monitor counts \n" "# \n" "# \n" "# \n" "#$$$$$$$.......... \n" " lambda 5.5e-10 0 0 0 x[ 0] \n" " P_polariz 1 0 0 0 x[ 1] \n" " N_polariz 0.5 0 0 0 x[ 2] \n" " P_analyse 1 0 0 0 x[ 3] \n" " N_analyse 0.5 0 0 0 x[ 4] \n" " amplitude 3.256 0.004329 0.1 1000 x[ 5] \n" " bg_const 3e-05 1e-07 0 1000 x[ 6] \n" " blur_i 0.01 0 0 0.2 x[ 7] \n" " blur_f 0.01 0 0 0.2 x[ 8] \n" " chsix 2e-07 0 0 0 x[ 9] \n" " chsiy 2e-07 0 0 0 x[10] \n" " chsiz 8e-05 0 0 0 x[11] \n" " B_enhance 1 0 0 0 x[12] \n" " G_enhance 1 0 0 0 x[13] \n" " alpha_foot 0.02 0 0 2 x[14] \n" " ai_offset 0 0 0 0 x[15] \n" " bg_factor 0 0 0 0 x[16] \n" " bg_alpha_i 0 0 0 0 x[17] \n" " bg_ampl 0 0 0 0 x[18] \n" " w2harmonic 0 0 0 0 x[19] \n" " absorber1 1 0 0 0 x[20] \n" " absorber2 1 0 0 0 x[21] \n" " absorber3 1 0 0 0 x[22] \n" " absorber4 1 0 0 0 x[23] \n" " absorber5 1 0 0 0 x[24] \n" " . 0 0 0 0 x[25] \n" " . 0 0 0 0 x[26] \n" " . 0 0 0 0 x[27] \n" " . 0 0 0 0 x[28] \n" " . 0 0 0 0 x[29] \n" " . 0 0 0 0 x[30] \n" " . 0 0 0 0 x[31] \n" " . 0 0 0 0 x[32] \n" " . 0 0 0 0 x[33] \n" " . 0 0 0 0 x[34] \n" " nb_vac_re 0 0 0 0 x[35] \n" " nb_vac_im 0 0 0 0 x[36] \n" " nb_bulk 1.32e+14 0 1.9e+13 2.1e+16 x[37] \n" " nb_imag 1e+10 0 0 2.1e+16 x[38] \n" " roughness 5e-10 0 0 5e-09 x[39] \n" "$ 40 4 #_group_repetition:_40_times;_this_is_the_No.4_of_4_in_the_group \n" " nb_real 5e+14 0 1.02e+14 2.8e+14 x[40] \n" " nb_imag 0 0 0 9e+14 x[41] \n" " d 2.587e-08 1.209e-10 2.5e-10 1e-06 x[42] \n" " roughness 0.1 0 0 0.25 x[43] \n" " B 0.8 0 0 0 x[44] \n" " fi 30 0 0 0 x[45] \n" " th 0 0 0 0 x[46] \n" "$ 40 3 #_group_repetition:_40_times;_this_is_the_No.3_of_4_in_the_group \n" " nb_real 9e+14 0 1.02e+14 2.8e+14 x[47] \n" " nb_imag 0 0 0 9e+14 x[48] \n" " d 2.587e-08 1.209e-10 2.5e-10 1e-06 x[49] \n" " roughness 0.12 0 0 0.25 x[50] \n" " B 0.95 0 0 0 x[51] \n" " fi -50 0 0 0 x[52] \n" " th 0 0 0 0 x[53] \n" "$ 40 2 #_group_repetition:_40_times;_this_is_the_No.2_of_4_in_the_group \n" " nb_real 3e+14 0 1.02e+14 2.8e+14 x[54] \n" " nb_imag 0 0 0 9e+14 x[55] \n" " d 2.587e-08 1.209e-10 2.5e-10 1e-06 x[56] \n" " roughness 0.08 0 0 0.25 x[57] \n" " B 1 0 0 0 x[58] \n" " fi 60 0 0 0 x[59] \n" " th 0 0 0 0 x[60] \n" "$ 40 1 #_group_repetition:_40_times;_this_is_the_No.1_of_4_in_the_group \n" " nb_real 1e+15 0 1.02e+14 2.8e+14 x[61] \n" " nb_imag 0 0 0 9e+14 x[62] \n" " d 2.587e-08 1.209e-10 2.5e-10 1e-06 x[63] \n" " roughness 0.15 0 0 0.25 x[64] \n" " B 0.7 0 0 0 x[65] \n" " fi -50 0 0 0 x[66] \n" " th 0 0 0 0 x[67] \n" "$ 1 1 #_group_repetition:___1_times;_this_is_the_No.1_of_1_in_the_group \n" " nb_real 2.934e+14 1.036e+10 1.02e+14 3.8e+14 x[68] \n" " nb_imag 0 0 0 9e+14 x[69] \n" " d 4.107e-08 5.142e-11 2.5e-10 1e-06 x[70] \n" " roughness 0.05 0 0 0.25 x[71] \n" " B 0 0 0 0 x[72] \n" " fi 0 0 0 0 x[73] \n" " th 0 0 0 0 x[74] \n" "#$$$$$$$.......... \n" "#this was the top (reflexion) layer next to vacuum \n" "# \n" ); fclose (file); file = fopen("sample.dat", "w"); fprintf(file, "# \n" "# This is the 'sample.dat' example file for the program 'superfit' \n" "# \n" "# Site for comments. \n" "# \n" "#motor reflected background monitor absorber \n" "#alpha_i intensity intensity intensity group \n" "#[degrees] [counts] [counts] [counts] \n" "0.0500 173392 1591 89024 1 \n" "0.1000 219335 2132 103048 1 \n" "0.2500 73186 736 145121 1 \n" "0.3000 28457 357 159146 1 \n" "0.1500 242465 2329 117073 1 \n" "0.2000 259944 2256 131097 1 \n" "0.3500 14173 278 173170 1 \n" "0.4000 7884 199 187195 1 \n" "0.4500 4750 188 201219 1 \n" "0.5000 3100 157 215243 1 \n" "0.5500 2116 163 229268 1 \n" "0.6000 1426 167 243292 2 \n" "0.6500 990 161 257317 2 \n" "0.7000 766 181 271341 2 \n" "0.7500 602 191 285365 2 \n" "0.8000 422 176 299390 2 \n" "0.8500 349 177 313414 2 \n" "0.9000 286 197 327439 2 \n" "0.9500 284 168 341463 2 \n" "1.0000 233 187 355487 2 \n" "1.0500 250 213 369512 2 \n" "1.1000 213 188 383536 2 \n" "1.1500 240 206 397560 2 \n" "1.2000 219 196 411585 2 \n" "1.2500 232 215 425609 2 \n" "1.3000 236 222 439634 2 \n" "1.3500 245 236 453658 2 \n" "1.4000 228 275 467682 2 \n" "1.4500 240 263 481707 2 \n" "1.5000 226 232 495731 2 \n" "1.5500 274 237 509756 2 \n" "1.6000 262 285 523780 2 \n" "1.6500 253 279 537804 2 \n" "1.7000 265 282 551829 2 \n" "1.7500 289 298 565853 2 \n" "1.8000 261 272 579878 2 \n" "1.8500 288 318 593902 2 \n" "1.9000 281 303 607926 2 \n" "1.9500 289 296 621951 2 \n" "2.0000 288 299 635975 2 \n" ); fclose (file); message ( "\nA minuit command sample file 'sample.min' has been written or rewritten.", "\nA data sample file 'sample.dat' has been written or rewritten."); message ("For a demonstration of the program in the 1D simulation mode, please type \n", " superfit sample 11 500 0 2 (it will take seconds). "); message ("For a demonstration of the program in the 1D fit mode, please type \n", " superfit sample 3 500 0.3 1.5 (it will take ~1 hour). "); message ("For a demonstration of the program in the 2D simulation mode, please type \n", " superfit sample 5 50 0 2 (it will take ~4 minutes). "); message ("For a demonstration of the program in the 2D simulation mode, please type \n", " superfit sample 5 500 0 5 (it will take hours). "); exit(1);} /********************** create sample files ******************/ /* Error functions: */ /*********************** terminate ****************/ void terminate(char *text){putchar('\a'); fflush(stdout); printf("%s\n",text); exit(1);} /************************* terminate ********************/ /********************** message ***************/ void message (char *text1, char *text2){ /*Text beim Suchen*/ /*putchar('\a'); fflush(stdout);*/ fprintf(stderr, "%s %s \n", text1, text2);} /********************** message ***************/ /********************** error ***************/ void error (char *ref){ /* reports on runtime errors usage: error ("label"); necessary: #include #include int errorcount = 0; int ERROR = 0; or 1 */ #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 "); message (text, " "); whistle(); if (ERROR) terminate (" Floating point error detected.");} feclearexcept (FE_ALL_EXCEPT); #endif /* fenv.h*/ errno = 0;} /********************** error ***************/ /********************** wait ********************/ void wait (char *text){ char ch[MAX_LEN]; putchar('\a'); fflush(stdout); fprintf(stderr, "Waiting for RETURN: %s", text); fgets(ch, MAX_LEN-1, stdin);} /********************** wait ********************/ /********************** whistle ********************/ void whistle (void){putchar('\a'); fflush(stdout);} /********************** whistle ********************/