diff --git a/README.md b/README.md index 275652d..a4b535b 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ Clara2 - a parallel classical radiation calculator based on Liénard-Wiechert po Introduction ------------ -Clara stands for CLAsical RAdiation and is a redevelopment from scratch of [Clara1](https://github.com/ComputationalRadiationPhysics/clara1). +Clara stands for CLAssical RAdiation and is a redevelopment from scratch of [Clara1](https://github.com/ComputationalRadiationPhysics/clara1). It has been developed as part of a [diploma thesis](http://www.hzdr.de/db/Cms?pOid=38997) in 2012. In contrast to [Clara1](https://github.com/ComputationalRadiationPhysics/clara1) it is parallelized using MPI and OpenMP to run efficiently on large CPU clusters. @@ -28,5 +28,14 @@ Software License ---------------- *Clara2* is licensed under the **GPLv3+**. You can use any of our *libraries* with -**GPLv3+ or LGPLv3+** (they are *dual licensed*). +**GPLv3+ or LGPLv3+** (they are *dual-licensed*). Please refer to our [LICENSE](LICENSE) + + +Dependency +---------- + +*Clara2* uses the FFTW library when used with the (faster) fft detector. +If you install *Clara2*, you need to install FFTW as well. You can finde +compiled code at: http://www.fftw.org/ or the source code at +https://github.com/FFTW/fftw3. FFTW3 is under GNU General Public License v2.0. diff --git a/REFERENCE.md b/REFERENCE.md index 844a41b..4a0ffad 100644 --- a/REFERENCE.md +++ b/REFERENCE.md @@ -10,3 +10,11 @@ Available online 13 November 2013, ISSN 0168-9002, http://dx.doi.org/10.1016/j.nima.2013.10.073. + +R. Pausch + +Electromagnetic Radiation from Relativistic Electrons as Characteristic Signature of their Dynamics + +Diploma Thesis 2012, Technische Universität Dresden, Germany + +DOI: 10.5281/zenodo.843510 \ No newline at end of file diff --git a/example.sh b/example.sh deleted file mode 100755 index be34ad4..0000000 --- a/example.sh +++ /dev/null @@ -1,20 +0,0 @@ -# Copyright 2014 Richard Pausch -# -# This file is part of Clara 2. -# -# Clara 2 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 3 of the License, or -# (at your option) any later version. -# -# Clara 2 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 should have received a copy of the GNU General Public License -# along with Clara 2. -# If not, see . -# - -echo "Hello World" diff --git a/src/all_directions.cpp b/src/all_directions.cpp index f0af90e..75df869 100644 --- a/src/all_directions.cpp +++ b/src/all_directions.cpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch, Joy, Alexander Koehler * * This file is part of Clara 2. * @@ -18,48 +18,27 @@ * If not, see . */ - - #include "all_directions.hpp" - - -#include -#include -#include -#include - -#include #include -#include // OpenMP - -#include "single_trace.hpp" - -#include "vector.hpp" -#include "analytical_solution.hpp" -#include "physics_units.hpp" -#include "string_manipulation.hpp" - -#include "import_from_file.hpp" +#include "single_direction.hpp" #include "load_txt.hpp" -#include "gzip_lib.hpp" - - - +#include "include/input_output.hpp" +#include "settings.hpp" +#include "setFilename.hpp" +#include "fileExists.hpp" /** * function that calculates spectra in different directions for * a single particle trace * - * @param trace_id a unique id which which the trajctopry file - * can be identified - * @param arg a string telling wether a "binary" or "asci" - * output should be used + * @param trace_id a unique id which identifies the trajectory file + * * @return error code **/ -int all_directions(const unsigned int trace_id, const char arg[]) +int all_directions(const unsigned int trace_id) { using namespace std; @@ -70,257 +49,187 @@ int all_directions(const unsigned int trace_id, const char arg[]) gettimeofday(&t1, NULL); - /* ------------ constants ------------------------------- */ - const double omega_max = 3.0e19; /* maximum of ploted frequency Hz */ - const double theta_max = 1.14594939; /* maximum of theta in degree */ - const unsigned int N_spectrum = 2048; /* number of frequencies "omega"*/ - const unsigned int N_theta = 120; /* number of directions in first angle "theta" */ - const unsigned int N_phi = 2; /* number of directions in second angle "phi" */ - const unsigned int N_trace = 2000; /* maximum number of traces */ - const unsigned int N_direction = N_theta*N_phi; // number of all directions + const unsigned int N_direction = param::N_theta * param::N_phi; // number of all directions /* ---------- get trace ID ----------------- */ /* check whether the id of the trace "trace_id" is larger than the given N_trace value */ - if(!(trace_id <= N_trace)) - { - std::cout << "trace-ID is out of range (MAX = " << N_trace << ")" << std::endl; - return 1; - } - - /* -------- get store info ----------- */ - /* since output can be stored as binary file or as asci file, here the selected option - * is checked or an error is thrown in case the selction was wrong */ - bool asci_output; - std::string store_str = arg; - if(!store_str.compare("asci")) - { - asci_output = true; - std::cout << "ASCI output" << std::endl; - } - else if(!store_str.compare("binary")) - { - asci_output = false; - std::cout << "binary output" << std::endl; - } - else - { - std::cerr << "2nd argument needs to be binary or asci" << std::endl; - throw "bin_asci"; - } - - - + if(!(trace_id <= param::N_trace)) + { + std::cout << "trace-ID is out of range (MAX = " << param::N_trace << ")" << std::endl; + return 1; + } /* ------- set up/compute all angles thetas ------------------ */ - double theta[N_theta]; - for(unsigned i=0; i< N_theta; ++i) - { - theta[i] = (double)i / N_theta * theta_max; - } + double theta[param::N_theta]; + for(unsigned i=0; i< param::N_theta; ++i) + { + theta[i] = (double)i / param::N_theta * param::theta_max; + } /* ------- set up/compute all angles phis ---------------------- */ - double phi[N_phi] = {0.0, 90.0}; - + double phi[param::N_phi] = {0.0, 90.0}; /* allocate memory for all spectra */ struct spectrum_container { - double spectrum[N_spectrum]; + double spectrum[param::N_spectrum]; }; - spectrum_container* all_spec = new spectrum_container[N_theta*N_phi]; - + spectrum_container* all_spec = new spectrum_container[N_direction]; /* compute the frequency array "omega" and fill spectra with zeros */ - double omega[N_spectrum]; - const double my_delta_omega = omega_max/N_spectrum; - for(unsigned i=0; i 254) - { - /* throw warning when buffer is to small for path name */ - std::cerr << "buffer to small!!! " << std::endl; - throw "Buffer to small!"; - } - /* print out path naem in order to check it in output files */ - std::cout << "check: filename: " << filename << std::endl; + char filenameTrace[param::N_char_filename]; + setFilename(filenameTrace, param::traceFileTemplate, trace_id, param::N_char_filename); + /* print out path name in order to check it in output files */ + std::cout << "check: filename: " << filenameTrace << std::endl; - /* -------- load trace from file ------- */ - + /* check if given file exists */ - if(!file_exists(filename)) + if(!file_exists(filenameTrace)) return 1; /* output to inform user that file is loaded */ - std::cout << "load file: " << filename << std::endl; + std::cout << "load file: " << filenameTrace << std::endl; /* create memory */ - const unsigned linenumber = linecounter(filename); /* get lines of data */ + const unsigned linenumber = linecounter(filenameTrace); /* get lines of data */ one_line* data = new one_line[linenumber]; /* get memory for data */ /* run function that fills data from file into "data": */ - load_txt(filename, linenumber, data); + load_txt(filenameTrace, linenumber, data); - - /* --------- calculate spectrum of one trace for all direction ---------- */ - /* in case of additional parallelsation using OpenMP uncomment this: */ + /* in case of additional parallelization using OpenMP remove comment: */ /* #pragma omp parallel for num_threads(4) schedule(dynamic, 1) */ for(unsigned direction_index = 0; direction_index< N_direction; ++direction_index) + { + const double my_theta = theta[direction_index % param::N_theta]; + const double my_phi = phi[direction_index/param::N_theta]; + printf("calculate direction: %4d -> theta: %3.5f , phi: %3.5f \n", direction_index, my_theta, my_phi); + + /* + * compute the spectra for a single direction + * and trow an error if something goes wrong + */ + if((single_direction(data, linenumber, omega, all_spec[direction_index].spectrum, param::N_spectrum, my_theta, my_phi))!=0) { - const double my_theta = theta[direction_index % N_theta]; - const double my_phi = phi[direction_index/N_theta]; - printf("calculate direction: %4d -> theta: %3.5f , phi: %3.5f \n", direction_index, my_theta, my_phi); - - /* - * compute the spectra for a single direction - * and trow an error if something goes wrong - */ - if((single_trace(data, linenumber, omega, all_spec[direction_index].spectrum, N_spectrum, my_theta, my_phi))!=0) - { - std::cerr << "error occured in single_trace function" << std::endl; - throw "error in single_trace function"; - } + std::cerr << "error occurred in single_direction function" << std::endl; + throw "error in single_direction function"; } - + } /* ------- outputfile ------------------------------ */ /* allocate memory for name of output file */ - char outputfilename[256]; - - /* fill output file for each trace_id based on template */ - if(sprintf(outputfilename, - "my_spectrum_trace%06d.dat", - trace_id) > 254) - { - std::cerr << "buffer to small!!! " << std::endl; - throw "Buffer to small!"; - } + char outputfilename[param::N_char_filename]; + setFilename(outputfilename, param::outputFileTemplate, trace_id, param::N_char_filename); /* print name of output file */ std::cout << "check: output-filename: " << outputfilename << std::endl; - /* --- file output ------------- */ - /* store spectral data either as binary or as asci data */ - if(asci_output) + /* store spectral data either as binary or as ascii data */ + if(param::ascii_output) + { + /* ---- ASCII output file ------------------------ */ + ofstream my_output(outputfilename); /* create file */ + if(my_output.is_open()) /* check if it is open */ { - /* ---- ASCI output file ------------------------ */ - ofstream my_output(outputfilename); /* create file */ - if(my_output.is_open()) /* check if it is open */ + for(unsigned j=0; j. */ - -#ifndef ALL_DIRECTIONS_RPAUSCH -#define ALL_DIRECTIONS_RPAUSCH +#pragma once /** * function that calculates spectra in different directions for * a single particle trace * - * @param trace_id a unique id which which the trajctopry file - * can be identified - * @param arg a string telling wether a "binary" or "asci" - * output should be used + * @param trace_id a unique id which which the trajectory file + * can be identified * @return error code **/ -int all_directions(const unsigned int trace_id, const char arg[]); - -#endif - +int all_directions(const unsigned int trace_id); diff --git a/src/clara2_hypnos.modules b/src/clara2_hypnos.modules index cc051a9..36f2145 100644 --- a/src/clara2_hypnos.modules +++ b/src/clara2_hypnos.modules @@ -1,5 +1,8 @@ - . /etc/profile.modules - module purge - module load gcc/4.6.2 - module load infiniband/1.0.0 - module load openmpi/1.6.0 +. /etc/profile.modules +module purge +module load gcc/4.6.2 +module load infiniband/1.0.0 +module load openmpi/1.6.0 +module load fftw/3.3.4 +module load editor/emacs + diff --git a/src/convert_to_matrix.cpp b/src/convert_to_matrix.cpp deleted file mode 100644 index 622c0f4..0000000 --- a/src/convert_to_matrix.cpp +++ /dev/null @@ -1,230 +0,0 @@ -/** - * Copyright 2014 Richard Pausch - * - * This file is part of Clara 2. - * - * Clara 2 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 3 of the License, or - * (at your option) any later version. - * - * Clara 2 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 should have received a copy of the GNU General Public License - * along with Clara 2. - * If not, see . - */ - - -/** - * This is a no longer used program's code that - * that was needed when a single MPI task calculated the - * radiation for only a single direction but for several - * particle traces. It collected all the results for - * different directions and created a single matrix-like - * output that contained the radiation spectra from all - * trajectories and all directions. - * - * SHOULD THIS STILL BE INCLUDED IN THE CODE? - ISSUE #17 - **/ - - - - -#include -#include -#include -#include -#include - - -/** TO DO: this should be in a separate file - ISSUE #15 **/ - -/** - * This function checks if a file exists on the hard drive. - * - * @param filename string containg the path and filename to bechecked - * @return Returns true if file exists, otherwise false. - **/ -bool file_exists(const char *filename) -{ - std::ifstream infile(filename); - return infile; -} - - - -/** - * This program combines files with spectra for differnt directions - * to a single files containing all spectra for directions. - * - * @return only 0 is return - no error code yet - **/ -int main() -{ - - /* verbose output to inform user that the program is running */ - std::cout << "Convert single file data to matrix: " << std::flush; - - /* ----------- set parameters ------------ */ - /* in this section, all parameters are set that need to be changed - * in order to set up the collection process. */ - - const unsigned N_omega = 2048; /* number of frequencies */ - /* To DO: could be read from linenumber input file or from a - * genral param file associated with the simulation - ISSUE #8 */ - - double omega[N_omega]; /* allocate memory for reading omega values in */ - - - - /* compute directions: theta and phi */ - /* The real values of theta and phi are needed because they are used - * for naming the files. */ - - - - /* compute directions: theta and phi */ - /* The real values of theta and phi are needed because they are used - * for naming the files. */ - - const unsigned N_theta = 120; /* number of different theta angles */ - /* To DO: this can not be read from input file but from a - * genral param file associated with the simulation - ISSUE #8 */ - - double theta[N_theta]; /* allocate memory for directions */ - - /* compute directions (theta): */ - double theta_max = 1.14594939; /* = 1.0/gamma in degree */ - for(unsigned i=0; i> omega[j] >> data_phi_0[i][j]; - } - else if(i> omega_temp >> data_phi_0[i][j]; - - /* check if previeously stored omegas are equal to - * those stored in the other files. If they are not - * equal write warning to screen. (but do not stop - * program) */ - if(omega_temp != omega[j]) - std::cout << " ERROR: frequency discrepancy: " << filename - << " - " << omega[j] << " != " << omega_temp << std::endl; - } - } - else /* for all files with phi=90 --> same procedure */ - { - double omega_temp; /* temporary variable for omega */ - /* go through data in files (see above's code */ - for(unsigned j=0; j> omega_temp >> data_phi_90[i%N_theta][j]; - if(omega_temp != omega[j]) - std::cout << " ERROR: frequency discrepancy: " << filename - << " - " << omega[j] << " != " << omega_temp << std::endl; - } - } - data.close(); /* close the file of thze current direction */ - } - } - - - /* -------- end load data --------- */ - - - /* store matrix-like spectral data for first phi - * tabs separate the values of differnt frequencies - * and newlines separate different directions. */ - std::ofstream output_0("matrix_phi_0.dat"); - for(unsigned i=0; i. - */ - - - -#include -#include -#include - -#ifndef GZIP_LIB_RPAUSCH -#define GZIP_LIB_RPAUSCH - - - -/* TO DO: not dealing with gzib files - ISSUE #19 */ -/** - * write any data as binary to file - * - * @param data pointer to data (array) - * @param size size of the data in bytes - * @param filename pointer to a char-array containing the filename - * to store data in - * @return gives zero on success, throws error if file can not - * be written to disk - **/ -int store_data(void* data, - long unsigned int size, - char* filename) -{ - /* create file handler to write in binary format = "wb" */ - FILE* outfile = fopen(filename, "wb"); - - if(outfile == NULL) /* in case file can not be written */ - { - std::cerr << "could not create file: " << filename << std::endl; - /* TO DO: error handling - ISSUE #20 */ - throw "could_not_create_file"; - } - - fwrite(data, 1, size, outfile); /* write data to file */ - fclose(outfile); /* close file handler */ - - return 0; -} - - - -/* TO DO: not dealing with gzib files - ISSUE #19 */ -/** - * read binary data from file - * - * @param data pointer to memory where data should be stored - * @param size number of bytes to be read from file and put into data - * @param filename pointer to char array with path and filename of - * the file containing the data - * @return zero on success, one in case the file can not be read - **/ -int read_data(void* data, - long unsigned int size, - char* filename) -{ - /* create file handler to read in binary format = "rb" */ - FILE* outfile = fopen(filename, "rb"); - - /* in case the file handler is erroneous (file does not exist, - * etc.) return 1 */ - if(outfile == NULL) - { - /* TO DO: error handling - ISSUE #20 */ - return 1; - } - - fread(data, 1, size, outfile); /* read data from file*/ - fclose(outfile); /* close file handler */ - - return 0; -} - - - -/* TO DO: is compresion realy used in the code? - ISSUE #5 */ -/** - * write data to file using compression - * - * @param data pointer to data that should be stored on disk - * @param size number of bytes that are in data and should be - * put in the file - * @param filename pointer to char-array containg path and - * filename of the file to put data in - * @param number "unsigned int" that is written to the - * beginning of the file - * @return returns zero in case of success, no error handling - **/ -int compress_data(void* data, - long unsigned int size, - char* filename, - unsigned int number) -{ - /* create compressed file handler (write binary="wb") */ - gzFile outfile = gzopen(filename, "wb"); - - /* write number in front of file */ - gzwrite(outfile, &number, sizeof(unsigned int)); - /* write data to file */ - gzwrite(outfile, data, size); - - /* close file handler */ - gzclose(outfile); - - return 0; -} - - - -/* TO DO: is compresion realy used in the code? - ISSUE #5 */ -/** - * append data to a compressed file - * - * @param data pointer to data that should be stored on disk - * @param size number of bytes that are in data and should be - * put in the file - * @param filename pointer to char-array containg path and - * filename of the file to append data to - * @return returns zero in case of success, no error handling - **/ -int compress_data_append(void* data, - long unsigned int size, - char* filename) -{ - /* create compressed file handler (append binary="ab") */ - gzFile outfile = gzopen(filename, "ab"); - - /* append data to file */ - gzwrite(outfile, data, size); - - /* close file handler */ - gzclose(outfile); - - return 0; -} - - - -/* TO DO: is this function still necessary? - ISSUE #5 */ -/** - * specific function to create a gzFile to access - * data from the "bubbleStressTest" simulation - * - * @param pFile empty "gzFile" that should point - * to simulation data - * @param od_or_even select whether the file lies in - * bigOutput1 (odd) or bigOutput2 - * (even) - * @param time_id simulation time step to identify - * file - * @param core_id CPU-core number for file name - **/ -void create_gzFile(gzFile& pFile, - int od_or_even, - unsigned int time_id, - unsigned int core_id) -{ - /* set up data file: */ - - /* create file name from "bubbleStressTest" path*/ - char filename[256]; - if(sprintf(filename, /* char-array to put path into */ - "/net/cns/projects/bubbleStressTest/bigOutput%1d/e%05d_%03d.dat.gz", - od_or_even, /* select path depending on odd or even */ - time_id, /* identify file by time step */ - core_id /* CPU core identification used for file naming */ - ) > 254) - { - /* throw error in case the buffer is to small */ - std::cerr << "buffer to small!!! " << std::endl; - /* TO DO: error handling - ISSUE #20 */ - throw "Buffer to small!"; - } - - /* create gzFile handler for read-only */ - pFile = gzopen(filename, "rb"); - - /* verify if file was created correctly */ - if(pFile == NULL) /* if error occured */ - { - pFile = gzopen(filename, "rb"); /* try again */ - if(pFile == NULL) // try again - { - /* if second attempt failes too throw error */ - std::cerr << "Could not open file" << std::endl; - /* TO DO: error handling - ISSUE #20 */ - throw "no file found"; - } - } -} - - - - - -#endif diff --git a/src/include/analytical_solution.hpp b/src/include/analytical_solution.hpp deleted file mode 100644 index 3bc43c5..0000000 --- a/src/include/analytical_solution.hpp +++ /dev/null @@ -1,97 +0,0 @@ -/** - * Copyright 2014 Richard Pausch - * - * This file is part of Clara 2. - * - * Clara 2 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 3 of the License, or - * (at your option) any later version. - * - * Clara 2 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 should have received a copy of the GNU General Public License - * along with Clara 2. - * If not, see . - */ - - - -#include -#include "utilities.hpp" -#include "physics_units.hpp" - -#ifndef ANALYTICAL_SOLUTION_RPAUSCH -#define ANALYTICAL_SOLUTION_RPAUSCH - - - -namespace analy -{ - - -const int N_per = 225; -const double lambda_u = 0.4E-6; -const double Bmax = 1.0; -const unsigned n_mode = 1; - -double K; -double F_approx; - - -using namespace std; - -double L(double x) -{ - return (util::square(sin(M_PI*x))) / (util::square(N_per*sin(M_PI*x/N_per))); - -} - -double L2(double x) -{ - return (util::square(sin(M_PI*x))) / (util::square(M_PI*x)); - -} - - - -double Lambda(double gamma) -{ - return lambda_u/(2*util::square(gamma)) * (1. + util::square(K)/2.); -} - -double Spectrum(double omega, double gamma) -{ - return util::square(phy::q*N_per*gamma)/ - (4.*M_PI*phy::epsilon_0*phy::c) * - L(N_per*(omega-2.*M_PI*phy::c/Lambda(gamma))/(2.*M_PI*phy::c / - Lambda(gamma)) ) - * F_approx; -} - -double Spectrum2(double omega, double gamma) -{ - return util::square(phy::q*N_per*gamma)/ - (4.*M_PI*phy::epsilon_0*phy::c) * - L2(N_per*(omega-2.*M_PI*phy::c/Lambda(gamma))/(2.*M_PI*phy::c / - Lambda(gamma)) ) - * F_approx; -} - - - -void calc_K(double gamma) -{ - K = 0.0569; //phy::q * Bmax * lambda_u / (2. * M_PI * phy::m_e * phy::c); - F_approx = util::square(K * n_mode) / util::square(1. + util::square(K)/2.) ; - - std::cout << " K= " << K << "\t F_approx = " << F_approx << std::endl; -} - - -} - -#endif diff --git a/src/include/detector_dft.cpp b/src/include/detector_dft.cpp index 0f79401..e963056 100644 --- a/src/include/detector_dft.cpp +++ b/src/include/detector_dft.cpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -22,155 +22,236 @@ #include "detector_dft.hpp" +#include "physics_units.hpp" +#include "utilities.hpp" + // Constructors and Destructors: -Detector_dft::Detector_dft(R_vec n_unit, double delta_t, const unsigned spek_length, - const double omega_max) - : n_unit(n_unit.unit_vec()), delta_t(delta_t), spek_length(spek_length), - spektrum(0), spektrum_mag(0), frequency(0) +/** constructor for a spectral detector + * + * @param n_unit = unit vector pointing in observation direction + * @param spec_length = number of frequencies to compute + * @param omega_max = maximum frequency -> frequency range [0, omega_max] + */ +Detector_dft::Detector_dft(const R_vec n_unit, + const unsigned spec_length, + const double omega_max) + : n_unit(n_unit.unit_vec()), + spec_length(spec_length), + spectrum(0), + spectrum_mag(0), + frequency(0) { - spektrum = new Vector, 3>[spek_length]; - spektrum_mag = new double[spek_length]; - frequency = new double[spek_length]; + spectrum = new Vector, 3>[spec_length]; /* complex amplitude for each omega */ + spectrum_mag = new double[spec_length]; /* absolute square of spectrum */ + frequency = new double[spec_length]; /* omega values */ - set_frequency(omega_max); + set_frequency(omega_max); /* compute frequencies */ } -Detector_dft::Detector_dft(R_vec n_unit, double delta_t, const unsigned spek_length, - const double* omega) - : n_unit(n_unit.unit_vec()), delta_t(delta_t), spek_length(spek_length), - spektrum(0), spektrum_mag(0), frequency(0) +/** constructor for a spectral detector + * + * @param n_unit = unit vector pointing in observation direction + * @param spec_length = number of frequenies to compute + * @param omega = pointer to frequencies to calculate + */ +Detector_dft::Detector_dft(const R_vec n_unit, + const unsigned spec_length, + const double* omega) + : n_unit(n_unit.unit_vec()), + spec_length(spec_length), + spectrum(0), spectrum_mag(0), + frequency(0) { - spektrum = new Vector, 3>[spek_length]; - spektrum_mag = new double[spek_length]; - frequency = new double[spek_length]; + spectrum = new Vector, 3>[spec_length]; /* complex amplitude for each omega */ + spectrum_mag = new double[spec_length]; /* absolute square of spectrum */ + frequency = new double[spec_length]; /* omega values */ - set_frequency(omega); + set_frequency(omega); /* copy frequency values */ } - +/** destructor */ Detector_dft::~Detector_dft() { - delete[] spektrum; - delete[] spektrum_mag; + delete[] spectrum; + delete[] spectrum_mag; delete[] frequency; } -// add to spectrum methodes: +// add to spectrum methods: -void Detector_dft::add_to_spectrum(const R_vec r, - const R_vec beta, - const R_vec dot_beta, - const double t_part) +/** method to add an radiation amplitude for one time step + * + * @param r_0 = position + * @param beta_0 = speed / speed_of_light + * @param dot_beta_0 = time derivative of beta_0 + * @param t_part_0 = time + * @param delta_t = time step width + */ +void Detector_dft::add_to_spectrum(const R_vec r, + const R_vec beta, + const R_vec dot_beta, + const double t_part, + const double delta_t) { - const R_vec fou1a = (n_unit%((n_unit-beta)%dot_beta)) - / util::square(1. - beta * n_unit); + /* calculate (real) vector part of radiation amplitude */ + const R_vec fou1 = (n_unit%((n_unit-beta)%dot_beta)) + / util::square(1. - beta * n_unit); - - const Vector, 3> fou1b = fou1a.make_complex(); + /* make vector part complex for calculation in complex space */ + const Vector, 3> fou1_complex = fou1.make_complex(); - for(unsigned a = 0; a < spek_length; ++a) + /* for each frequency: */ + for(unsigned a = 0; a < spec_length; ++a) + { + const double t_ret = t_part - (n_unit*r)/phy::c; /* retarded time */ + /* compute complex phase for this retarded time */ + std::complex fou2 = std::polar(1.0, frequency[a]*t_ret); + /* for each dimension: */ + for (unsigned i=0; i<3; ++i) { - std::complex fou2 = std::polar(1.0, frequency[a]* - (t_part - (n_unit*r)/phy::c)); - for (unsigned i=0; i<3; ++i) - { - (spektrum[a])[i] += fou1b[i]*fou2*delta_t; - } - + /* this is the (local) Nyquist limiter - only add amplitude if + * frequency is below (local) Nyquist frequency */ + if(frequency[a] < 0.75 * M_PI / (delta_t*(1.0 - beta * n_unit)) ) + (spectrum[a])[i] += fou1_complex[i]*fou2*delta_t; } + + } } -void Detector_dft::add_to_spectrum(const R_vec r, - const R_vec p, - const R_vec dot_p, - const R_vec beta, - const double gamma, - const double dot_gamma, - const double t_part) -{ +/* ISSUE # 74 - interface too crowded */ +/** method to add an radiation amplitude for one time step + * + * @param r_0 = position + * @param p_0 = momentum + * @param dot_p_0 = time derivative of p_0 + * @param beta_0 = speed / speed_of_light + * @param gamma_0 = relativistic gamma factor + * @param dot_gamma_0 = time derivative of gamma_0 + * @param t_part_0 = time + * @param delta_t = time step width + */ +void Detector_dft::add_to_spectrum(const R_vec r, + const R_vec p, + const R_vec dot_p, + const R_vec beta, + const double gamma, + const double dot_gamma, + const double t_part, + const double delta_t) +{ + /* normalize momentum */ const R_vec p_wave = p/(phy::m_e*phy::c); const R_vec dot_p_wave = dot_p/(phy::m_e*phy::c); + /* calculate (real) vector part of radiation amplitude */ const R_vec fou1 = (n_unit%((gamma*n_unit - p_wave)%(dot_p_wave - beta*dot_gamma))) - / util::square(gamma - p_wave*n_unit); + / util::square(gamma - p_wave*n_unit); + /* make vector part complex for calculation in complex space */ const Vector, 3> fou1_complex = fou1.make_complex(); - for(unsigned a = 0; a < spek_length; ++a) + /* for each frequency: */ + for(unsigned a = 0; a < spec_length; ++a) + { + const double t_ret = t_part - (n_unit*r)/phy::c; /* retarded time */ + /* compute complex phase for this retarded time */ + std::complex fou2 = std::polar(1.0, frequency[a]*t_ret); + /* for each dimension: */ + for (unsigned i=0; i<3; ++i) { - std::complex fou2 = std::polar(1.0, frequency[a]* - (t_part - (n_unit*r)/phy::c)); - for (unsigned i=0; i<3; ++i) - { - (spektrum[a])[i] += fou1_complex[i]*fou2*delta_t; - } - + /* this is the (local) Nyquist limiter - only add amplitude if + * frequency is below (local) Nyquist frequency */ + if(frequency[a] < 0.75 * M_PI / (delta_t*(1.0 - beta * n_unit)) ) + (spectrum[a])[i] += fou1_complex[i]*fou2*delta_t; } + + } } -// calculate spectrum methode: +// calculate spectrum method: +/** calculate the spectrum from all added amplitudes */ void Detector_dft::calc_spectrum() { - const double factor = util::square(phy::q)/ - (16.*util::cube(M_PI)*phy::epsilon_0*phy::c); - - for(unsigned a = 0; a < spek_length; ++a) - { - spektrum_mag[a] = factor * - util::square(spektrum[a].abs()); - } + /* SI unit factor for calculating radiation energy per frequency and solid angle */ + const double factor = util::square(phy::q) + / (16.*util::cube(M_PI)*phy::epsilon_0*phy::c); + + /* calculate the absolute value for each frequency */ + for(unsigned a = 0; a < spec_length; ++a) + { + spectrum_mag[a] = factor * util::square(spectrum[a].abs()); + } } // Getter: -double Detector_dft::get_spectrum(unsigned a, unsigned b) +/* ISSUE #75 - separate get spectrum and get frequency */ +/** get spectrum or frequency + * + * @param a = id between [0, number of frequencies] + * @param b = int return 0->frequency, 1->spectra in Js + */ +double Detector_dft::get_spectrum(unsigned a, + unsigned b) { - assert(a < spek_length); + assert(a < spec_length); switch(b) - { - case 0: + { + case 0: /* return frequency */ return frequency[a]; break; - case 1: - return spektrum_mag[a]; + case 1: /* return spectral value at that frequency */ + return spectrum_mag[a]; break; - default: - std::cout << "Wrong access to spektrum (dft, beta)!" << std::endl; + default: /* index out of range */ + std::cout << "Wrong access to spectrum (dft, beta)!" << std::endl; std::cout << b << " is larger than 2." << std::endl; assert(false); break; - } + } } +/** return total energy (in calculated spectral range) */ double Detector_dft::energy() { - double result = 0.; - for (unsigned i = 0; i < (spek_length); ++i) - result += spektrum_mag[i]; + double result = 0.; /* total energy */ + /* sum over all frequencies */ + for (unsigned i = 0; i < (spec_length); ++i) + result += spectrum_mag[i]; + + /* multiply by frequency bin width */ + /* ISSUE #76 - remove this arbitrary delta_omega here */ result *= (frequency[7] - frequency[6]); return result; } -// set frequency methodes: +// set frequency methods: -inline void Detector_dft::set_frequency(const double omega_max) +/** set frequency based on maximum frequency + * to cover frequency range [0, omega_max] + * + * @ param omega_max = maximum frequency + */ +inline void Detector_dft::set_frequency(const double omega_max) { - for (unsigned i=0; i -#include -#include -#include +#pragma once #include "vector.hpp" -#include "physics_units.hpp" -#include "utilities.hpp" - -#ifndef DETECTOR_DFT_RPAUSCH -#define DETECTOR_DFT_RPAUSCH - -//! \brief class for a point-like detector storing the signal externaly +/** Detector class that computes the spectra via a + * (non-equidistat) discrete Fourier transform + */ class Detector_dft { public: - //! \brief constructor for a spectral detector - /*! - @param n_unit = unit vector in direction of energy deposition - @param delta_t = timestep of odint - */ - Detector_dft(R_vec n_unit, double delta_t, const unsigned spek_length, - const double omega_max); - - Detector_dft(R_vec n_unit, double delta_t, const unsigned spek_length, - const double* omega); - - - ~Detector_dft(); - - void add_to_spectrum(const R_vec r_0, - const R_vec beta_0, - const R_vec dot_beta_0, - const double t_part_0); - - void add_to_spectrum(const R_vec r_0, - const R_vec p_0, - const R_vec dot_p_0, - const R_vec beta_0, - const double gamma_0, - const double dot_gamma_0, - const double t_part_0); - + /** constructor for a spectral detector + * + * @param n_unit = unit vector pointing in observation direction + * @param spec_length = number of frequenies to compute + * @param omega_max = maximum frequency -> frequency range [0, omega_max] + */ + Detector_dft(const R_vec n_unit, + const unsigned spec_length, + const double omega_max); + + /** constructor for a spectral detector + * + * @param n_unit = unit vector pointing in observation direction + * @param spec_length = number of frequenies to compute + * @param omega = pointer to frequencies to calculate + */ + Detector_dft(const R_vec n_unit, + const unsigned spec_length, + const double* omega); + + + /** destructor */ + ~Detector_dft(); + + /** method to add an radiation amplitude for one time step + * + * @param r_0 = position + * @param beta_0 = speed / speed_of_light + * @param dot_beta_0 = time derivative of beta_0 + * @param t_part_0 = time + * @param delta_t = time step width + */ + void add_to_spectrum(const R_vec r_0, + const R_vec beta_0, + const R_vec dot_beta_0, + const double t_part_0, + const double delta_t); + + /* ISSUE # 74 - interface too crowded */ + /** method to add an radiation amplitude for one time step + * + * @param r_0 = position + * @param p_0 = momentum + * @param dot_p_0 = time derivative of p_0 + * @param beta_0 = speed / speed_of_light + * @param gamma_0 = relativistic gamma factor + * @param dot_gamma_0 = time derivative of gamma_0 + * @param t_part_0 = time + * @param delta_t = time step width + */ + void add_to_spectrum(const R_vec r_0, + const R_vec p_0, + const R_vec dot_p_0, + const R_vec beta_0, + const double gamma_0, + const double dot_gamma_0, + const double t_part_0, + const double delta_t); + + /** calculate the spectrum from all added amplitudes */ void calc_spectrum(); - double get_spectrum(unsigned a, unsigned b); - - double energy(); - + /** get spectrum or frequency + * + * @param a = id between [0, number of frequencies] + * @param b = int return 0->frequency, 1->spectra in Js + */ + double get_spectrum(unsigned a, + unsigned b); -private: - // data: - const R_vec n_unit; - const double delta_t; - const unsigned spek_length; - + /** return total energy (in calculated spectral range) */ + double energy(); -public: // debugging - Vector, 3>* spektrum; private: - double* spektrum_mag; - double* frequency; - - - // set frequency methodes: + // data: + const R_vec n_unit; /* observation direction */ + const unsigned spec_length; /* number of frequencies to calculate */ + Vector, 3>* spectrum; /* complex amplitudes */ + double* spectrum_mag; /* resulting spectrum */ + double* frequency; /* pointer to frequency values */ + + // set frequency methods: + + /** set frequency based on maximum frequency + * to cover frequency range [0, omega_max] + * + * @ param omega_max = maximum frequency + */ inline void set_frequency(const double omega_max); + + /** set frequency using pointer to frequency values + * + * @ param omega = double pointer to frequency values + */ inline void set_frequency(const double* omega); }; - - -#endif - diff --git a/src/include/detector_e_field.cpp b/src/include/detector_e_field.cpp index ef6aeb9..0875a5c 100644 --- a/src/include/detector_e_field.cpp +++ b/src/include/detector_e_field.cpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,39 +19,54 @@ */ +#include "detector_e_field.hpp" +#include "physics_units.hpp" +#include "utilities.hpp" -#include "detector_e_field.hpp" -//! \brief calculates retardated signal at detector -/*! @param r_0 = position one time step before last step - @param r_1 = position at last time step - @param p_0 = momentum one time step before last step - @param p_1 = momentum at last time step - @param dot_p_0 = dp/dt one time step before last step - @param dot_p_1 = dp/dt at last time step - @param beta_0 = v/c = beta one time step befor last step - @param beta_1 = v/c = beta at last time step - @param gamma_0 = gamma one time step befor last step - @param gamma_1 = gamma at last time step - @param dot_gamma_0 = d gamma/dt one time step befor last step - @param dot_gamma_1 = d gamma/dt at last time step - @param t_part_0 = time at the electron one time step before last step - */ - -inline Detector_e_field::Detector_e_field(R_vec detector, double delta_t, unsigned N_sig, - unsigned start_sig ) - : delta_t(delta_t), signal(N_sig, start_sig), detector(detector), - counter(0) +/** constructor for electric field detector (at position in space) + * + * @param detector = R_vec with observation position + * @param delta_t = time step width + * @param N_sig = maximum number of field entries to the detector + * @param start_sig = start time at which the electric field should + * be "recorded" + */ +inline Detector_e_field::Detector_e_field(R_vec detector, + double delta_t, + unsigned N_sig, + unsigned start_sig) + : delta_t(delta_t), + signal(N_sig, start_sig), + detector(detector), + counter(0) { } - + +/** compute integer index in which the retarded electric field (signal) + * should be written to + * + * @param t_signal = time at which the electric field arrives at + * detector position + * @return unsigned int index in which the time fits best in + * the signal class + */ inline unsigned Detector_e_field::delay_index(double t_signal) { return unsigned(t_signal / delta_t) +1 ; } - -inline double Detector_e_field::t_signal(R_vec r, double t) + +/** compute the time at which the electric field emitted by + * a charged particle at position r and time t arrives at the + * detector position (retarded time) + * + * @param r = R_vec position of particle + * @param t = double current time (of the particle) + * @return signal arrival time (retarded time) + */ +inline double Detector_e_field::t_signal(R_vec r, + double t) { R_vec verbindung = r - detector; // (? do I use this vector ?) double R = verbindung.mag(); @@ -59,147 +74,198 @@ inline double Detector_e_field::t_signal(R_vec r, double t) } - +/** returns the counter of double counts minus no counts + * + * @return value of index counter + */ inline int Detector_e_field::count() -{ - return counter; +{ + return counter; } +/* ISSUE #74 - clean up interface */ +/** compute electric field value(s) at detector position + * and store the field values using two time steps + * (subscript 0 and 1) + * + * @param r_0 = particle position at t_part_0 + * @param r_1 = particle position at t_part_0 + delta_t + * @param p_0 = particle momentum at t_part_0 + * @param p_1 = particle momentum at t_part_0 + delta_t + * @param dot_p_0 = particle change in momentum at t_part_0 + * @param dot_p_1 = particle change in momentum at t_part_0 + delta_t + * @param beta_0 = particle beta (v/c) at t_part_0 + * @param beta_1 = particle beta (v/c) at t_part_0 + delta_t + * @param gamma_0 = particle relativistic gamma at t_part_0 + * @param gamma_1 = particle relativistic gamma at t_part_0 + delta_t + * @param dot_gamma_0 = particle change in gamma at t_part_0 + * @param dot_gamma_1 = particle change in gamma at t_part_0 + delta_t + * @param t_part_0 = time at 'first' time step + */ +void Detector_e_field::place(const R_vec r_0, + const R_vec r_1, + const R_vec p_0, + const R_vec p_1, + const R_vec dot_p_0, + const R_vec dot_p_1, + const R_vec beta_0, + const R_vec beta_1, + const double gamma_0, + const double gamma_1, + const double dot_gamma_0, + const double dot_gamma_1, + const double t_part_0) +{ + /* compute electric field (signal) arrival times at detector + * and map this to a discrete index entry in the detector electric + * field array */ + double t_signal_0 = t_signal(r_0, t_part_0); + double t_signal_1 = t_signal(r_1, t_part_0+delta_t); + unsigned delay_index_0 = delay_index(t_signal_0); + unsigned delay_index_1 = delay_index(t_signal_1); -void Detector_e_field::place(const R_vec r_0, const R_vec r_1, - const R_vec p_0, const R_vec p_1, - const R_vec dot_p_0, const R_vec dot_p_1, - const R_vec beta_0, const R_vec beta_1, - const double gamma_0, const double gamma_1, - const double dot_gamma_0, const double dot_gamma_1, - const double t_part_0) -{ - // signal arrival times discrete - double t_signal_0 = t_signal(r_0, t_part_0); - double t_signal_1 = t_signal(r_1, t_part_0+delta_t); - - unsigned delay_index_0 = delay_index(t_signal_0); - unsigned delay_index_1 = delay_index(t_signal_1); - - if (delay_index_0 == delay_index_1 -1) - { - // preparation for Lienard Wiechert - double t_prim = delta_t * delay_index_0 - t_signal_0; - R_vec r_prim = interpol(r_0, r_1, t_prim); - R_vec p_prim = interpol(p_0, p_1, t_prim); - R_vec dot_p_prim = interpol(dot_p_0, dot_p_1, t_prim); - - double gamma_prim = interpol(gamma_0, gamma_1, t_prim); - double dot_gamma_prim = interpol(dot_gamma_0, dot_gamma_1, t_prim); - R_vec beta_prim = interpol(beta_0, beta_1, t_prim); - R_vec beta_dot_times_gamma_prim = (dot_p_prim*(1./(phy::m_e*phy::c)) - - dot_gamma_prim*beta_prim); - - R_vec verbindung_prim = r_prim - detector; - R_vec e_R_prim = verbindung_prim.unit_vec(); - double R_prim = verbindung_prim.mag(); - - - signal[delay_index_0] = Lienard_Wiechert(e_R_prim, p_prim, - beta_dot_times_gamma_prim, - gamma_prim, R_prim); - } - - else if (delay_index_0 == delay_index_1) - { - counter-- ; - } - - else if (delay_index_0 +2 == delay_index_1 ) - { - counter++ ; - // preparation for Lienard Wiechert - double t_prim1 = delta_t * delay_index_0 - t_signal_0; - R_vec r_prim1 = interpol(r_0, r_1, t_prim1); - R_vec p_prim1 = interpol(p_0, p_1, t_prim1); - R_vec dot_p_prim1 = interpol(dot_p_0, dot_p_1, t_prim1); - - double gamma_prim1 = interpol(gamma_0, gamma_1, t_prim1); - double dot_gamma_prim1 = interpol(dot_gamma_0, dot_gamma_1, t_prim1); - R_vec beta_prim1 = interpol(beta_0, beta_1, t_prim1); - R_vec beta_dot_times_gamma_prim1 = (dot_p_prim1*(1./(phy::m_e*phy::c)) - - dot_gamma_prim1*beta_prim1); - - R_vec verbindung_prim1 = r_prim1 - detector; - R_vec e_R_prim1 = verbindung_prim1.unit_vec(); - double R_prim1 = verbindung_prim1.mag(); - - // --------------------------- - - double t_prim2 = delta_t * (delay_index_0+1) - t_signal_0; - R_vec r_prim2 = interpol(r_0, r_1, t_prim2); - R_vec p_prim2 = interpol(p_0, p_1, t_prim2); - R_vec dot_p_prim2 = interpol(dot_p_0, dot_p_1, t_prim2); - - double gamma_prim2 = interpol(gamma_0, gamma_1, t_prim2); - double dot_gamma_prim2 = interpol(dot_gamma_0, dot_gamma_1, t_prim2); - R_vec beta_prim2 = interpol(beta_0, beta_1, t_prim2); - R_vec beta_dot_times_gamma_prim2 = (dot_p_prim2*(1./(phy::m_e*phy::c)) - - dot_gamma_prim2*beta_prim2); - - R_vec verbindung_prim2 = r_prim2 - detector; - R_vec e_R_prim2 = verbindung_prim2.unit_vec(); - double R_prim2 = verbindung_prim2.mag(); - - - signal[delay_index_0] = Lienard_Wiechert(e_R_prim1, p_prim1, - beta_dot_times_gamma_prim1, - gamma_prim1, R_prim1); - signal[delay_index_0+1] = Lienard_Wiechert(e_R_prim2, p_prim2, - beta_dot_times_gamma_prim2, - gamma_prim2, R_prim2); - } - - else - { - std::cout << "Unbekannter Fehler " << std::endl; - assert(false); - } -} + /* ISSUE #78 - time index might be off by one */ + /* in case both time step would be closest to two neighboring + * array entries interpolate between both time step to fill one + * array entry */ + if (delay_index_0 == delay_index_1 -1) + { + /* time difference between discrete index and retarded time */ + double t_prim = delta_t * delay_index_0 - t_signal_0; + /* interpolate between time step 0 and 1 */ + R_vec r_prim = interpol(r_0, r_1, t_prim); + R_vec p_prim = interpol(p_0, p_1, t_prim); + R_vec dot_p_prim = interpol(dot_p_0, dot_p_1, t_prim); + double gamma_prim = interpol(gamma_0, gamma_1, t_prim); + double dot_gamma_prim = interpol(dot_gamma_0, dot_gamma_1, t_prim); + R_vec beta_prim = interpol(beta_0, beta_1, t_prim); + /* calculate beta * gamma at time t_prim (+ t_signal_0) */ + R_vec beta_dot_times_gamma_prim = (dot_p_prim*(1./(phy::m_e*phy::c)) + - dot_gamma_prim*beta_prim); + /* calculate distance between particle and detector position and derive + * necessary quantities */ + R_vec distance_prim = r_prim - detector; + R_vec e_R_prim = distance_prim.unit_vec(); + double R_prim = distance_prim.mag(); + /* ISSUE #78 - time index might be off by one */ + /* write electric field into electric field array (signal) */ + signal[delay_index_0] = Lienard_Wiechert(e_R_prim, p_prim, + beta_dot_times_gamma_prim, + gamma_prim, R_prim); + } -//! \brief simple interpolation between two values (f.e. location, speed) -/*! @param r_0 = value at startpoint - @param r_1 = value at endpoint (one timestep later) - @param t = time as interpolation parameter (0 < t < delta_t) - */ -template -V Detector_e_field::interpol(V r_0, V r_1, double t) -// interpolation between 2 points -{ - return r_0 + ((r_1 - r_0)*(t/delta_t)); -} + /* both time steps would be closest to the same time index + * do not calculate the electric field */ + else if (delay_index_0 == delay_index_1) + { + counter--; /* reduce counter by one */ + } + + /* ISSUE #78 - time index might be off by one */ + /* the two time step cover two three discrete time steps + * do multiple interpolations */ + else if (delay_index_0 +2 == delay_index_1 ) + { + counter++ ; + /* time difference between discrete index and retarded time */ + double t_prim1 = delta_t * delay_index_0 - t_signal_0; + /* ---- first time index ---- */ + /* interpolate between time step 0 and 1 for first index */ + R_vec r_prim1 = interpol(r_0, r_1, t_prim1); + R_vec p_prim1 = interpol(p_0, p_1, t_prim1); + R_vec dot_p_prim1 = interpol(dot_p_0, dot_p_1, t_prim1); + double gamma_prim1 = interpol(gamma_0, gamma_1, t_prim1); + double dot_gamma_prim1 = interpol(dot_gamma_0, dot_gamma_1, t_prim1); + R_vec beta_prim1 = interpol(beta_0, beta_1, t_prim1); + /* calculate beta * gamma at time t_prim1 (+ t_signal_0) */ + R_vec beta_dot_times_gamma_prim1 = (dot_p_prim1*(1./(phy::m_e*phy::c)) + - dot_gamma_prim1*beta_prim1); + /* calculate distance between particle and detector position and derive + * necessary quantities */ + R_vec distance_prim1 = r_prim1 - detector; + R_vec e_R_prim1 = distance_prim1.unit_vec(); + double R_prim1 = distance_prim1.mag(); + /* ISSUE #78 - time index might be off by one */ + /* write electric field into electric field array (signal) */ + signal[delay_index_0] = Lienard_Wiechert(e_R_prim1, p_prim1, + beta_dot_times_gamma_prim1, + gamma_prim1, R_prim1); + /* ---- second time index ---- */ + /* interpolate between time step 0 and 1 for second index */ + double t_prim2 = delta_t * (delay_index_0+1) - t_signal_0; + R_vec r_prim2 = interpol(r_0, r_1, t_prim2); + R_vec p_prim2 = interpol(p_0, p_1, t_prim2); + R_vec dot_p_prim2 = interpol(dot_p_0, dot_p_1, t_prim2); + double gamma_prim2 = interpol(gamma_0, gamma_1, t_prim2); + double dot_gamma_prim2 = interpol(dot_gamma_0, dot_gamma_1, t_prim2); + R_vec beta_prim2 = interpol(beta_0, beta_1, t_prim2); + /* calculate beta * gamma at time t_prim2 (+ t_signal_0) */ + R_vec beta_dot_times_gamma_prim2 = (dot_p_prim2*(1./(phy::m_e*phy::c)) + - dot_gamma_prim2*beta_prim2); + /* calculate distance between particle and detector position and derive + * necessary quantities */ + R_vec distance_prim2 = r_prim2 - detector; + R_vec e_R_prim2 = distance_prim2.unit_vec(); + double R_prim2 = distance_prim2.mag(); + + /* ISSUE #78 - time index might be off by one */ + /* write electric field into electric field array (signal) */ + signal[delay_index_0+1] = Lienard_Wiechert(e_R_prim2, p_prim2, + beta_dot_times_gamma_prim2, + gamma_prim2, R_prim2); + } + /* in case none of the above occurs - throw an error + * (if more discrete time steps are in between the linear interpolation + * would hide the electric field dynamic due to a bad temporal resolution) */ + else + { + std::cout << "unknown error " << std::endl; + assert(false); + } +} -/*! \brief calculates Lienard Wiechert Potation (more precise the - $\vec E$-field */ -/*! @param e_R = unit vector in the from the electron to the detector - @param beta = beta vector of the electron - $ \vec \beta = \frac{\vec v}{c} $ - @param beta_dot = $ \operatorname{\frac{d}{dt}} \vec \beta $ - @param gamma = gamma value of the electron - $ gamma = \sqrt{\frac{1}{1-\vec \beta^2} } $ - @ param R = distance between electron and detector - */ -R_vec Detector_e_field::Lienard_Wiechert(const R_vec& e_R, const R_vec& p, - const R_vec beta_dot_times_gamma, const double& gamma, - const double& R) -{ - R_vec p_o = p*(1/(phy::c*phy::m_e)); - - return (1./(4.*M_PI*phy::epsilon_0))* - ( phy::q * e_R % ((gamma*e_R-p_o) % beta_dot_times_gamma)*gamma / - (phy::c*R*util::cube(gamma-p_o*e_R)) ) ; +/** simple interpolation between two values (f.e. location, speed) + * at two time step with time difference delta_t + * + * @param x_0 = value at start point + * @param x_1 = value at endpoint (one time step later) + * @param t = time as interpolation parameter (0 < t < delta_t) + */ +template +V Detector_e_field::interpol(V x_0, + V x_1, + double t) +// interpolation between 2 points +{ + return x_0 + ((x_1 - x_0)*(t/delta_t)); } +/** calculates the electric field $\vec E$ based on Lienard Wiechert potential + * @param e_R = unit vector pointing from the electron to the detector + * @param beta = beta of the electron + * $ \vec \beta = \frac{\vec v}{c} $ + * @param beta_dot = $ \operatorname{\frac{d}{dt}} \vec \beta $ + * @param gamma = gamma value of the electron + * $ gamma = \sqrt{\frac{1}{1-\vec \beta^2} } $ + * @param R = distance between electron and detector + * @return electric field at detector position + */ +R_vec Detector_e_field::Lienard_Wiechert(const R_vec& e_R, + const R_vec& p, + const R_vec beta_dot_times_gamma, + const double& gamma, + const double& R) +{ + R_vec p_o = p*(1/(phy::c*phy::m_e)); + return (1./(4.*M_PI*phy::epsilon_0)) + * ( phy::q * e_R % ((gamma*e_R-p_o) % beta_dot_times_gamma)*gamma + / (phy::c*R*util::cube(gamma-p_o*e_R)) ) ; +} diff --git a/src/include/detector_e_field.hpp b/src/include/detector_e_field.hpp index 706cf81..0b2e492 100644 --- a/src/include/detector_e_field.hpp +++ b/src/include/detector_e_field.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,106 +19,127 @@ */ - - -#include -#include -#include -#include +#pragma once #include "vector.hpp" #include "large_index_storage.hpp" -#include "physics_units.hpp" -#include "utilities.hpp" - - - -#ifndef DETECTOR_E_FIELD_RPAUSCH -#define DETECTOR_E_FIELD_RPAUSCH -//! \brief class for a point-like detector storing the signal externaly +/** \brief class for a point-like detector storing the electric field (signal) externally */ class Detector_e_field { public: - //! \brief constructor for a point-like detector - /*! @param detector = location of the detector - @param delta_t = timestep of odint - @param N_sig = number of datapoints of signal to store - @param start_sig = start index of signal - signal = pointer to signal struct array (signal at detector) - */ - inline Detector_e_field(R_vec detector, double delta_t, unsigned N_sig, - unsigned start_sig ); - - - //! \brief Calculate delay_index from signal time - //! @param t_signal Time at which the signal arrives at the detector - inline unsigned delay_index(double t_signal); - - - //! \brief Calculate time of signal at the detector - /*! @param r position of moving charge (electron) - @param t time at which the particle moved */ - inline double t_signal(R_vec r, double t); - - - // for details see .cpp file: - void place(const R_vec r_0, const R_vec r_1, - const R_vec p_0, const R_vec p_1, - const R_vec dot_p_0, const R_vec dot_p_1, - const R_vec beta_0, const R_vec beta_1, - const double gamma_0, const double gamma_1, - const double dot_gamma_0, const double dot_gamma_1, - const double t_part_0); - - - /// returns the counter of double counts minus no counts + /** constructor for electric field detector (at position in space) + * + * @param detector = R_vec with observation position + * @param delta_t = time step width + * @param N_sig = maximum number of field entries to the detector + * @param start_sig = start time at which the electric field should + * be "recorded" + */ + inline Detector_e_field(R_vec detector, + double delta_t, + unsigned N_sig, + unsigned start_sig ); + + /** compute integer index in which the retarded electric field (signal) + * should be written to + * + * @param t_signal = time at which the electric field arrives at + * detector position + * @return unsigned int index in which the time fits best in + * the signal class + */ + inline unsigned delay_index(double t_signal); + + /** compute the time at which the electric field emitted by + * a charged particle at position r and time t arrives at the + * detector position (retarded time) + * + * @param r = R_vec position of particle + * @param t = double current time (of the particle) + * @return signal arrival time (retarded time) + */ + inline double t_signal(R_vec r, double t); + + /* ISSUE #74 - clean up interface */ + /** compute electric field value(s) at detector position + * and store the field values using two time steps + * (subscript 0 and 1) + * + * @param r_0 = particle position at t_part_0 + * @param r_1 = particle position at t_part_0 + delta_t + * @param p_0 = particle momentum at t_part_0 + * @param p_1 = particle momentum at t_part_0 + delta_t + * @param dot_p_0 = particle change in momentum at t_part_0 + * @param dot_p_1 = particle change in momentum at t_part_0 + delta_t + * @param beta_0 = particle beta (v/c) at t_part_0 + * @param beta_1 = particle beta (v/c) at t_part_0 + delta_t + * @param gamma_0 = particle relativistic gamma at t_part_0 + * @param gamma_1 = particle relativistic gamma at t_part_0 + delta_t + * @param dot_gamma_0 = particle change in gamma at t_part_0 + * @param dot_gamma_1 = particle change in gamma at t_part_0 + delta_t + * @param t_part_0 = time at 'first' time step + */ + void place(const R_vec r_0, + const R_vec r_1, + const R_vec p_0, + const R_vec p_1, + const R_vec dot_p_0, + const R_vec dot_p_1, + const R_vec beta_0, + const R_vec beta_1, + const double gamma_0, + const double gamma_1, + const double dot_gamma_0, + const double dot_gamma_1, + const double t_part_0); + + /** returns the counter of double counts minus no counts + * + * @return value of index counter + */ inline int count(); - -// data: - const double delta_t; /// length of timesteps - Large_index_storage signal; /// E_field at detector - - -private: - - //data - R_vec detector; /// location of the detector - int counter; - - - //methodes - - //! \brief simple interpolation between two values (f.e. location, speed) - /*! @param r_0 = value at startpoint - @param r_1 = value at endpoint (one timestep later) - @param t = time as interpolation parameter (0 < t < delta_t) - */ - template - V interpol(V r_0, V r_1, double t); // interpolation between - // 2 points - - /*! \brief calculates Lienard Wiechert Potation (more precise the - $\vec E$-field */ - /*! @param e_R = unit vector in the from the electron to the detector - @param beta = beta vector of the electron - $ \vec \beta = \frac{\vec v}{c} $ - @param beta_dot = $ \operatorname{\frac{d}{dt}} \vec \beta $ - @param gamma = gamma value of the electron - $ gamma = \sqrt{\frac{1}{1-\vec \beta^2} } $ - @ param R = distance between electron and detector - */ - R_vec Lienard_Wiechert(const R_vec& e_R, const R_vec& p, - const R_vec beta_dot_times_gamma, - const double& gamma, - const double& R); - - - -}; + /* data: */ + const double delta_t; /* length of time steps */ + Large_index_storage signal; /* container for electric field at detector */ -#endif +private: + /* data: */ + R_vec detector; /* location of the detector */ + int counter; /* the counter of double counts minus no counts + * (info for user to optimze detector) not needed for calulations */ + + /* methods: */ + + /** simple interpolation between two values (f.e. location, speed) + * at two time step with time difference delta_t + * + * @param x_0 = value at start point + * @param x_1 = value at endpoint (one time step later) + * @param t = time as interpolation parameter (0 < t < delta_t) + */ + template + V interpol(V r_0, + V r_1, + double t); + + /** calculates the electric field $\vec E$ based on Lienard Wiechert potential + * @param e_R = unit vector pointing from the electron to the detector + * @param beta = beta of the electron + * $ \vec \beta = \frac{\vec v}{c} $ + * @param beta_dot = $ \operatorname{\frac{d}{dt}} \vec \beta $ + * @param gamma = gamma value of the electron + * $ gamma = \sqrt{\frac{1}{1-\vec \beta^2} } $ + * @param R = distance between electron and detector + * @return electric field at detector position + */ + R_vec Lienard_Wiechert(const R_vec& e_R, + const R_vec& p, + const R_vec beta_dot_times_gamma, + const double& gamma, + const double& R); +}; diff --git a/src/include/detector_fft.cpp b/src/include/detector_fft.cpp index 197b34a..025711b 100644 --- a/src/include/detector_fft.cpp +++ b/src/include/detector_fft.cpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch, Alexander Koehler * * This file is part of Clara 2. * @@ -19,138 +19,193 @@ */ - - #include "detector_fft.hpp" +#include "../settings.hpp" +#include "physics_units.hpp" +#include "utilities.hpp" +#include "ned_fft.hpp" + // Constructor and Destructor: -Detector_fft::Detector_fft(R_vec n_unit, unsigned N_data) - : n_unit(n_unit.unit_vec()), N_data(N_data), counter(0), - time(0), data(0), spektrum_mag(0), frequency(0) + +/** constructor for a point-like detector using a Fast Fourier Transform + * to calculate the radiation spectra + * + * @param n_unit = unit vector in direction of energy deposition + * @param N_data = number of time step in trace + */ +Detector_fft::Detector_fft(const R_vec n_unit, + const unsigned N_data) + : n_unit(n_unit.unit_vec()), + N_data(N_data), + counter(0), + time(0), + data(0), + spectrum_mag(0), + frequency(0) { - spek_length = power_of_two(N_data); + spec_length = power_of_two(N_data) * param::fft_length_factor; - //std::cout << "spek_length : " << spek_length << std::endl; - time = new double[spek_length]; - data = new R_vec[spek_length]; - spektrum_mag = new double[spek_length]; - frequency = new double[spek_length]; + time = new double[spec_length]; /* retarded time */ + data = new R_vec[spec_length]; /* real vector amplitude */ + spectrum_mag = new double[spec_length]; /* absolute square of spectrum */ + frequency = new double[spec_length]; /* omega values */ } +/** destructor */ Detector_fft::~Detector_fft() { delete[] time; delete[] data; - delete[] spektrum_mag; + delete[] spectrum_mag; delete[] frequency; } +// add to spectrum methods: + +/** method to add an radiation amplitude for one time step + * + * @param r = position + * @param beta = speed / speed_of_light + * @param dot_beta = time derivative of beta + * @param t_part = time + * @param delta_t = time step width + */ +void Detector_fft::add_to_spectrum(const R_vec r, + const R_vec beta, + const R_vec dot_beta, + const double t_part, + const double delta_t) +{ + /* calculate (real) vector part of radiation amplitude */ + const R_vec fou1 = (n_unit%((n_unit-beta)%dot_beta)) + / util::cube(1. - beta * n_unit); + /* abort if trace is longer than spectra */ + assert(counter < spec_length); // --> still necessary? probably not -// add to spectrum methodes: - -void Detector_fft::add_to_spectrum(const R_vec r, - const R_vec beta, - const R_vec dot_beta, - const double t_part, - const double delta_t) -{ - const R_vec fou1 = (n_unit%((n_unit-beta)%dot_beta)) - / util::cube(1. - beta * n_unit); - - assert(counter < spek_length); // --> still necessary? probably not - + /* retarded time */ time[counter] = t_part - (n_unit*r)/phy::c; + /* real vector amplitude */ data[counter] = fou1; - ++counter; + ++counter; /* count processed time steps */ } -void Detector_fft::add_to_spectrum(const R_vec r, - const R_vec p, - const R_vec dot_p, - const R_vec beta, - const double gamma, - const double dot_gamma, - const double t_part, - const double delta_t) -{ - +/* ISSUE # 74 - interface too crowded */ +/** method to add an radiation amplitude for one time step + * + * @param r = position + * @param p = momentum + * @param dot_p = time derivative of p + * @param beta = speed / speed_of_light + * @param gamma = relativistic gamma factor + * @param dot_gamma = time derivative of gamma + * @param t_part = time + * @param delta_t = time step width + */ +void Detector_fft::add_to_spectrum(const R_vec r, + const R_vec p, + const R_vec dot_p, + const R_vec beta, + const double gamma, + const double dot_gamma, + const double t_part, + const double delta_t) +{ + /* normalize momentum */ const R_vec p_wave = p/(phy::m_e*phy::c); const R_vec dot_p_wave = dot_p/(phy::m_e*phy::c); - const R_vec fou = (n_unit%((gamma*n_unit - p_wave)%(dot_p_wave - - beta*dot_gamma))) / util::cube(gamma - p_wave*n_unit) - * gamma; - - assert(counter < spek_length); // --> still necessairy? probably not + /* calculate (real) vector part of radiation amplitude */ + const R_vec fou = (n_unit%((gamma*n_unit - p_wave)%(dot_p_wave + - beta*dot_gamma))) / util::cube(gamma - p_wave*n_unit) + * gamma; + /* abort if trace is longer than spectra */ + assert(counter < spec_length); // --> still necessary? probably not + + /* retarded time */ time[counter] = t_part - (n_unit*r)/phy::c; + /* real vector amplitude */ data[counter] = fou; - ++counter; + ++counter; /* count processed time steps */ } -// calculate spectrum methode: +// calculate spectrum method: +/** calculate the spectrum from all added amplitudes */ void Detector_fft::calc_spectrum() { - for (unsigned i= counter; i analyse(spek_length, time, data); - //std::cout << " delta_t_ret = " << analyse.delta_t << std::endl; + /* SI unit factor for calculating radiation energy per frequency and solid angle */ + const double factor = util::square(phy::q) + / (16.*util::cube(M_PI)*phy::epsilon_0*phy::c); + /* compute spectra using a not-equidistant FFT method */ + ned_FFT analyse(spec_length, time, data); - for(unsigned i=0; ifrequency, 1->spectra in Js + */ double Detector_fft::get_spectrum(unsigned a, unsigned b) { - assert(a < (spek_length) ); + assert(a < (spec_length) ); switch(b){ - case 0: + case 0: /* return frequency */ return frequency[a]; - break; //neccesairy? - case 1: - return spektrum_mag[a]; + break; //necessary? + case 1: /* return spectral value at that frequency */ + return spectrum_mag[a]; break; - default: - std::cerr << "Wrong access to spektrum (fft, momentum)!" << std::endl; + default: /* index out of range */ + std::cerr << "Wrong access to spectrum (fft, momentum)!" << std::endl; std::cerr << b << " is larger than 1." << std::endl; assert(false); break; } } - +/** return total energy (in calculated spectral range) */ double Detector_fft::energy() { - double result = 0; + double result = 0; /* total energy */ + /* sum over all frequencies */ for (unsigned i = 0; i < half_frequency(); ++i) - result += spektrum_mag[i]; + result += spectrum_mag[i]; + + /* multiply by frequency bin width */ + /* ISSUE #76 - remove this arbitrary delta_omega here */ result *= (frequency[7] - frequency[6]); return result; } +/** return the number of frequencies till the Nyquist frequency + * (number of half of all frequency bins) + */ unsigned Detector_fft::half_frequency() { - return spek_length>>1; + return spec_length>>1; /* = spec_length / 2 */ } - diff --git a/src/include/detector_fft.hpp b/src/include/detector_fft.hpp index bd0f81e..91988c2 100644 --- a/src/include/detector_fft.hpp +++ b/src/include/detector_fft.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,84 +19,97 @@ */ - - -#include -#include -#include -#include +#pragma once #include "vector.hpp" -#include "physics_units.hpp" -#include "utilities.hpp" -#include "fft_ned.hpp" - -#ifndef DETECTOR_FFT_RPAUSCH -#define DETECTOR_FFT_RPAUSCH - -//! \brief class for a point-like detector storing the signal externaly +/** \brief class for a point-like detector storing the signal externally */ class Detector_fft { public: - //! \brief constructor for a point-like detector - /*! - @param n_unit = unit vector in direction of energy deposition - @param delta_t = timestep of odint - */ - Detector_fft(R_vec n_unit, unsigned N_data); - + /** constructor for a point-like detector using a Fast Fourier Transform + * to calculate the radiation spectra + * + * @param n_unit = unit vector in direction of energy deposition + * @param N_data = number of time step in trace + */ + Detector_fft(const R_vec n_unit, + const unsigned N_data); + + /** destructor */ ~Detector_fft(); - void add_to_spectrum(const R_vec r_0, - const R_vec beta_0, - const R_vec dot_beta_0, - const double t_part_0, - const double delta_t); - - - void add_to_spectrum(const R_vec r_0, - const R_vec p_0, - const R_vec dot_p_0, - const R_vec beta_0, - const double gamma_0, - const double dot_gamma_0, - const double t_part_0, - const double delta_t); - - + /** method to add an radiation amplitude for one time step + * + * @param r = position + * @param beta = speed / speed_of_light + * @param dot_beta = time derivative of beta + * @param t_part = time + * @param delta_t = time step width + */ + void add_to_spectrum(const R_vec r_0, + const R_vec beta_0, + const R_vec dot_beta_0, + const double t_part_0, + const double delta_t); + + + /* ISSUE # 74 - interface too crowded */ + /** method to add an radiation amplitude for one time step + * + * @param r = position + * @param p = momentum + * @param dot_p = time derivative of p + * @param beta = speed / speed_of_light + * @param gamma = relativistic gamma factor + * @param dot_gamma = time derivative of gamma + * @param t_part = time + * @param delta_t = time step width + */ + void add_to_spectrum(const R_vec r_0, + const R_vec p_0, + const R_vec dot_p_0, + const R_vec beta_0, + const double gamma_0, + const double dot_gamma_0, + const double t_part_0, + const double delta_t); + + /** calculate the spectrum from all added amplitudes */ void calc_spectrum(); - + /* ISSUE #75 - separate get spectrum and get frequency */ + /** get spectrum or frequency + * + * @param a = id between [0, number of frequencies] + * @param b = int return 0->frequency, 1->spectra in Js + */ double get_spectrum(unsigned a, unsigned b); + /** return total energy (in calculated spectral range) */ double energy(); + /** return the number of frequencies till the Nyquist frequency + * (number of half of all frequency bins) + */ unsigned half_frequency(); +//data: private: - //data: - - const R_vec n_unit; - // delta_t; // neccesairy for integration - const unsigned N_data; + const R_vec n_unit; /* observation direction */ + const unsigned N_data; /* number of time steps in trace */ - //public: // --> better to private !!! - unsigned spek_length; - unsigned counter; + unsigned spec_length; /* number of frequencies for spectra */ + unsigned counter; /* number of time steps analyzed */ - double* time; - R_vec* data; + double* time; /* pointer to retarded time */ + R_vec* data; /* pointer to real vector amplitudes */ - double* spektrum_mag; + double* spectrum_mag; /* pointer to absolute spectra */ public: - double* frequency; - -}; - - -#endif + double* frequency; /* pointer to frequencies */ +}; diff --git a/src/include/discrete.hpp b/src/include/discrete.hpp index 1968682..4719b8c 100644 --- a/src/include/discrete.hpp +++ b/src/include/discrete.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,161 +19,206 @@ */ - - +#pragma once #include "utilities.hpp" #include "physics_units.hpp" -#ifndef DISCRETE_RPAUSCH -#define DISCRETE_RPAUSCH - - /** - * \brief: storage class to handle 4 values and calculate derivatives - * usage Descrete - */ - + * \brief: storage class to handle 4 values and calculate derivatives + * usage Discrete + */ template class Discrete { - // friend class Retardation; public: - - //! \brief constructor filling all values - /*! @param old2 value at time = t-3 - @param old value at time = t-2 - @param now value at time = t-1 - @param future value at time = t - @param h pointer to discrete of time */ - Discrete(T old2, T old, T now, T future, const Discrete* h = 0) - : old2(old2), old(old), now(now), future(future), h(h) {} - - //! \brief constructor without filling data - /*! @param h difference between time steps */ - Discrete(const Discrete* h) - : h(h) {} - - //! \brief copy constructor - Discrete& operator=(const Discrete& copy) - { - assert(copy.h == h); - old2 = copy.old2; - old = copy.old; - now = copy.now; - future = copy.future; - - return *this; - } - - //! \brief setting new futre value and moving data down (now --> old, ...) - //! @param next new future value - void next(T next) - { - old2 = old; - old = now; - now = future; - future = next; - } - - //! \brief returning derivative for time = t - 2 - T dot_old() const - {return (now - old2)/(h->get_now() - h->get_old2()); } - - //! \brief returning derivative for time = t - 1 - T dot_now() const - {return (future - old)/(h->get_future() - h->get_old()); } - - //! \brief return value for time = t - 3 (old2) - T get_old2() const - {return old2; } - - //! \brief return value for time = t - 2 (old) - T get_old() const - {return old; } - - //! \brief return value for time = t - 1 (now) - T get_now() const - {return now; } - - //! \brief return value for time = t - 0 (future) - T get_future() const - {return future; } - - //! \brief return value for time = t - 0 (future) - T get_delta_old() const - {return (get_now()-get_old()); } - - + + /** \brief constructor filling all values + * + * @param old2 value at time = t-3 + * @param old value at time = t-2 + * @param now value at time = t-1 + * @param future value at time = t + * @param h pointer to discrete time + */ + Discrete(T old2, + T old, + T now, + T future, + const Discrete* h = 0) + : old2(old2), + old(old), + now(now), + future(future), + h(h) + { } + + /** \brief constructor without filling data + * + * @param h difference between time steps + */ + Discrete(const Discrete* h = 0) + : h(h) + { } + + /** \brief copy constructor + * + * @param copy Discrete object to copy data from + * @return the object itself + */ + Discrete& operator=(const Discrete& copy) + { + assert(copy.h == h); + old2 = copy.old2; + old = copy.old; + now = copy.now; + future = copy.future; + + return *this; + } + + /** \brief setting new future value and moving data down (now --> old, ...) + * + * @param next new future value + */ + void next(T next) + { + old2 = old; + old = now; + now = future; + future = next; + } + + /** \brief returning derivative for time = t - 2 + */ + T dot_old() const + { + /* second order symmetric time derivative */ + return (now - old2)/(h->get_now() - h->get_old2()); + } + + /** \brief returning derivative for time = t - 1 + */ + T dot_now() const + { + /* second order symmetric time derivative */ + return (future - old)/(h->get_future() - h->get_old()); + } + + /** \brief return value for time = t - 3 (old2) + */ + T get_old2() const + { + return old2; + } + + /** \brief return value for time = t - 2 (old) + */ + T get_old() const + { + return old; + } + + /** \brief return value for time = t - 1 (now) + */ + T get_now() const + { + return now; + } + + /** \brief return value for time = t - 0 (future) + */ + T get_future() const + { + return future; + } + + /** \brief return value for time = t - 0 (future) + */ + T get_delta_old() const + { + return (get_now()-get_old()); + } + + private: - T old2, old, now, future; - const Discrete* h; + T old2; /* value at t-3 */ + T old; /* value at t-2 */ + T now; /* value at t-1 */ + T future; /* value at t-0 */ + const Discrete* h; /* pointer to time values */ }; -/*! - * \brief a class provinding additional methods for relativistic physics\n - * returns Discrete objects for gamma and beta\n - * and also provides functors for calculating single values of gamma and beta\n - * Energie = sqrt(p^2*c^2+m_0^2*c^4) = GAMMA*m - * \vec beta = (\vec v)/c with \vec v = SPEED = (\vec p)/m(v) - * = (\vec p)/(m_0 * gamma) - */ +/** \brief a class providing additional methods for relativistic physics\n + * returns Discrete objects for gamma and beta\n + * and also provides functors for calculating single values of gamma and beta\n + * Energy = sqrt(p^2*c^2+m_0^2*c^4) = GAMMA*m + * \vec beta = (\vec v)/c with \vec v = SPEED = (\vec p)/m(v) + * = (\vec p)/(m_0 * gamma) + */ class More_discrete { public: - //! \brief constructor - /*! @param Det is a reference to a detector class from which one gets: - -> stepwidth: the length of the timestep between to steps */ - More_discrete(Discrete* h) - : stepwidth(h) {} - - - //! \brief returns a Discrete object for gamma - /*! @param p is a Discrete object of the momentum */ - Discrete momentum_to_gamma(Discrete& p) - { - return Discrete(gamma(p.get_old2()), - gamma(p.get_old()), - gamma(p.get_now()), - gamma(p.get_future()), - stepwidth); - } - - //! \brief returns a Discrete object for beta - /*! @param p is a Discrete object of the momentum - @param gamma is a Discrete object of gamma --> better to calculate here? - */ - Discrete momentum_to_beta(Discrete p, Discrete gamma) - { - return Discrete(beta(p.get_old2(), gamma.get_old2()), - beta(p.get_old(), gamma.get_old()), - beta(p.get_now(), gamma.get_now()), - beta(p.get_future(), gamma.get_future()), - stepwidth); - } - - //! \brief functor to calculate a single gamma value - //! @param p momentum - double gamma(R_vec p){ - return sqrt(util::square(p*phy::c)+ - util::square(phy::m_e*util::square(phy::c)))/(phy::m_e* - util::square(phy::c)); - } - - //! \brief functor to calculate a single beta value - /*! @param p momentum - @param gamma gamma */ - R_vec beta(R_vec p, double gamma) - { - return p*(1.0/(phy::c*phy::m_e*gamma)); - } - -private: - const Discrete* stepwidth; - -}; + /** \brief constructor + * + * @param h is a pointer to the time values + */ + More_discrete(const Discrete* h) + : stepwidth(h) {} + + + /** \brief convert momentum to gamma values + * + * @param p is a Discrete object of the momentum + */ + Discrete momentum_to_gamma(Discrete& p) + { + return Discrete(gamma(p.get_old2()), + gamma(p.get_old()), + gamma(p.get_now()), + gamma(p.get_future()), + stepwidth); + } + + /** \brief converts momentum to relativistic beta (v/c) + * + * @param p is a Discrete object of the momentum + * @param gamma is a Discrete object of gamma --> better to calculate here? + */ + Discrete momentum_to_beta(Discrete p, + Discrete gamma) + { + return Discrete(beta(p.get_old2(), gamma.get_old2()), + beta(p.get_old(), gamma.get_old()), + beta(p.get_now(), gamma.get_now()), + beta(p.get_future(), gamma.get_future()), + stepwidth); + } + + /** \brief method to calculate a single gamma value from a given momentum value + * + * @param p momentum + */ + double gamma(R_vec p) + { + return sqrt(util::square(p*phy::c)+ + util::square(phy::m_e*util::square(phy::c))) + /(phy::m_e*util::square(phy::c)); + } + + /** \brief method to calculate a single beta value from a given momentum and gamma value + * + * @param p momentum + * @param gamma relativistic gamma + */ + R_vec beta(R_vec p, + double gamma) + { + return p*(1.0/(phy::c*phy::m_e*gamma)); + } +private: + const Discrete* stepwidth; -#endif - +}; diff --git a/src/include/fileExists.cpp b/src/include/fileExists.cpp new file mode 100644 index 0000000..a22cfdc --- /dev/null +++ b/src/include/fileExists.cpp @@ -0,0 +1,35 @@ +/** + * Copyright 2016 Richard Pausch + * + * This file is part of Clara 2. + * + * Clara 2 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 3 of the License, or + * (at your option) any later version. + * + * Clara 2 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 should have received a copy of the GNU General Public License + * along with Clara 2. + * If not, see . + */ + +#include "fileExists.hpp" + +#include + + +/** \brief check whether a file exists or not + * + * @param filename pointer to array containing file location + * @return Returs true if file exists, otherwise false. + */ +bool file_exists(const char *filename) +{ + std::ifstream ifile(filename); + return ifile; +} diff --git a/example.cpp b/src/include/fileExists.hpp similarity index 72% rename from example.cpp rename to src/include/fileExists.hpp index 95900aa..d2464b5 100644 --- a/example.cpp +++ b/src/include/fileExists.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2016 Richard Pausch * * This file is part of Clara 2. * @@ -18,7 +18,13 @@ * If not, see . */ -int main() -{ - return 0; -} + +#pragma once + + +/** \brief check whether a file exists or not + * + * @param filename pointer to array containing file location + * @return Returs true if file exists, otherwise false. + */ +bool file_exists(const char *filename); diff --git a/src/include/import_from_file.hpp b/src/include/import_from_file.hpp index a05a1ee..440c6fd 100644 --- a/src/include/import_from_file.hpp +++ b/src/include/import_from_file.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,19 +19,14 @@ */ +#pragma once -#include - -#ifndef IMPORT_FROM_FILE_RPAUSCH -#define IMPORT_FROM_FILE_RPAUSCH - - -//! \brief simple container to store data from the Clara trace -/*! usage: one_line x[number of data lines]; then x[i].intern_data[0-6] */ +/* ISSUE #84 - re-factor this struct */ +/** \brief simple container to store data from the loaded traces + * usage: one_line x[number of data lines]; then x[i].intern_data[0-6] + * e.g. x[time_step].intern_data[0] = position_x; + */ struct one_line { - double intern_data[7]; // simple data structur two handel 7 doubles per line + double intern_data[7]; /* simple data structure to handle 7 doubles per line */ }; - -#endif - diff --git a/src/include/input_output.hpp b/src/include/input_output.hpp new file mode 100644 index 0000000..f21712e --- /dev/null +++ b/src/include/input_output.hpp @@ -0,0 +1,87 @@ +/** + * Copyright 2014-2016 Richard Pausch + * + * This file is part of Clara 2. + * + * Clara 2 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 3 of the License, or + * (at your option) any later version. + * + * Clara 2 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 should have received a copy of the GNU General Public License + * along with Clara 2. + * If not, see . + */ + + +#pragma once + +#include + + +/** + * write any data as binary to file + * + * @param data pointer to data (array) + * @param size size of the data in bytes + * @param filename pointer to a char-array containing the filename + * to store data in + * @return gives zero on success, throws error if file can not + * be written to disk + **/ +int store_data(void* data, + long unsigned int size, + char* filename) +{ + /* create file handler to write in binary format = "wb" */ + FILE* outfile = fopen(filename, "wb"); + + if(outfile == NULL) /* in case file can not be written */ + { + std::cerr << "could not create file: " << filename << std::endl; + /* TO DO: error handling - ISSUE #20 */ + throw "could_not_create_file"; + } + + fwrite(data, 1, size, outfile); /* write data to file */ + fclose(outfile); /* close file handler */ + + return 0; +} + + + +/** + * read binary data from file + * + * @param data pointer to memory where data should be stored + * @param size number of bytes to be read from file and put into data + * @param filename pointer to char array with path and filename of + * the file containing the data + * @return zero on success, one in case the file can not be read + **/ +int read_data(void* data, + long unsigned int size, + char* filename) +{ + /* create file handler to read in binary format = "rb" */ + FILE* outfile = fopen(filename, "rb"); + + /* in case the file handler is erroneous (file does not exist, + * etc.) return 1 */ + if(outfile == NULL) + { + /* TO DO: error handling - ISSUE #20 */ + return 1; + } + + fread(data, 1, size, outfile); /* read data from file*/ + fclose(outfile); /* close file handler */ + + return 0; +} diff --git a/src/include/interpolation.hpp b/src/include/interpolation.hpp index 88e20ce..0d161de 100644 --- a/src/include/interpolation.hpp +++ b/src/include/interpolation.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,32 +19,38 @@ */ +#pragma once -#ifndef INTERPOLATION_RPAUSCH -#define INTERPOLATION_RPAUSCH - -#include "detector_fft.hpp" template -void interpolation(const X* x_old, const Y* y_old, const unsigned N_old, - const X* x_new, Y* y_new, const unsigned N_new); - +void interpolation(const X* x_old, + const Y* y_old, + const unsigned N_old, + const X* x_new, + Y* y_new, + const unsigned N_new); template -void interpolation_on(const X* x_old, const Y* y_old, const unsigned N_old, - const X* x_new, Y* y_new, const unsigned N_new); +void interpolation_on(const X* x_old, + const Y* y_old, + const unsigned N_old, + const X* x_new, + Y* y_new, + const unsigned N_new); -template -void interpolation_on(const Detector_fft* fft, - const X* x_new, Y* y_new, const unsigned N_new); +template +void interpolation_on(const Detector_fft* fft, + const X* x_new, + Y* y_new, + const unsigned N_new); -void interpolation_int(Detector_fft* fft, - const double* x_new, double* y_new, const unsigned N_new); +void interpolation_int(Detector_fft* fft, + const double* x_new, + double* y_new, + const unsigned N_new); #include "interpolation.tpp" - -#endif diff --git a/src/include/interpolation.tpp b/src/include/interpolation.tpp index 9a129f0..87dbbdd 100644 --- a/src/include/interpolation.tpp +++ b/src/include/interpolation.tpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,117 +19,134 @@ */ -#include "interpolation.hpp" template -void interpolation(const X* x_old, const Y* y_old, const unsigned N_old, - const X* x_new, Y* y_new, const unsigned N_new) +void interpolation(const X* x_old, + const Y* y_old, + const unsigned N_old, + const X* x_new, + Y* y_new, + const unsigned N_new) { for(unsigned i_new = 0, i_old = 0; i_new< N_new; ++i_new) + { + while(x_old[i_old] < x_new[i_new] && i_old < N_old) + ++i_old; + if(i_old < N_old && i_old > 0) { - while(x_old[i_old] < x_new[i_new] && i_old < N_old) - ++i_old; - if(i_old < N_old && i_old > 0) - { - y_new[i_new] = (Y)((y_old[i_old] - y_old[i_old -1])/(x_old[i_old] - x_old[i_old -1]) * - (x_new[i_new] - x_old[i_old -1]) + y_old[i_old -1]); - } - else - y_new[i_new] = 0; - } + y_new[i_new] = (Y)((y_old[i_old] - y_old[i_old -1])/(x_old[i_old] - x_old[i_old -1]) + * (x_new[i_new] - x_old[i_old -1]) + y_old[i_old -1]); + } + else + { + y_new[i_new] = 0; + } + } } - - template -void interpolation_on(const X* x_old, const Y* y_old, const unsigned N_old, - const X* x_new, Y* y_new, const unsigned N_new) +void interpolation_on(const X* x_old, + const Y* y_old, + const unsigned N_old, + const X* x_new, + Y* y_new, + const unsigned N_new) { for(unsigned i_new = 0, i_old = 0; i_new< N_new; ++i_new) + { + while(x_old[i_old] < x_new[i_new] && i_old < N_old) + ++i_old; + if(i_old < N_old && i_old > 0) + { + y_new[i_new] += (Y)((y_old[i_old] - y_old[i_old -1])/(x_old[i_old] - x_old[i_old -1]) + * (x_new[i_new] - x_old[i_old -1]) + y_old[i_old -1]); + } + else { - while(x_old[i_old] < x_new[i_new] && i_old < N_old) - ++i_old; - if(i_old < N_old && i_old > 0) - { - y_new[i_new] += (Y)((y_old[i_old] - y_old[i_old -1])/(x_old[i_old] - x_old[i_old -1]) * - (x_new[i_new] - x_old[i_old -1]) + y_old[i_old -1]); - } - else - y_new[i_new] += (Y)0; // could be skipped - } + y_new[i_new] += (Y)0; // could be skipped + } + } } - -void interpolation_on(Detector_fft* fft, - const double* x_new, double* y_new, const unsigned N_new) +void interpolation_on(Detector_fft* fft, + const double* x_new, + double* y_new, + const unsigned N_new) { unsigned N_old = fft->half_frequency(); for(unsigned i_new = 0, i_old = 0; i_new< N_new; ++i_new) + { + while(fft->get_spectrum(i_old, 0) < x_new[i_new] && i_old < N_old) + ++i_old; + if(i_old < N_old && i_old > 0) + { + y_new[i_new] += ((fft->get_spectrum(i_old, 1) - fft->get_spectrum(i_old-1, 1))/(fft->get_spectrum(i_old, 0) - fft->get_spectrum(i_old-1, 0)) + * (x_new[i_new] - fft->get_spectrum(i_old-1, 0)) + fft->get_spectrum(i_old-1, 1)); + } + else { - while(fft->get_spectrum(i_old, 0) < x_new[i_new] && i_old < N_old) - ++i_old; - if(i_old < N_old && i_old > 0) - { - y_new[i_new] += ((fft->get_spectrum(i_old, 1) - fft->get_spectrum(i_old-1, 1))/(fft->get_spectrum(i_old, 0) - fft->get_spectrum(i_old-1, 0)) * - (x_new[i_new] - fft->get_spectrum(i_old-1, 0)) + fft->get_spectrum(i_old-1, 1)); - } - else - y_new[i_new] += 0; // could be skipped - } + y_new[i_new] += 0; // could be skipped + } + } } -void interpolation_int(Detector_fft* fft, - const double* x_new, double* y_new, const unsigned N_new) + +void interpolation_int(Detector_fft* fft, + const double* x_new, + double* y_new, + const unsigned N_new) { unsigned N_old = fft->half_frequency(); unsigned total_counter =0; unsigned j=0; - //std::cout << " started integrated interpolation: " << N_old << std::endl; - //in between: for(unsigned i=0; iget_spectrum(j, 0) >= x_low && fft->get_spectrum(j, 0) <= x_high && jget_spectrum(j, 1); - ++j; - } - - y_new[i] += sum/counter; - //std::cout << " ->" << counter << std::flush; - total_counter += counter; + x_high = 0.5 * (x_new[i] + x_new[i+1]); + x_low = x_new[i]; + } + else if(i==N_new-1) + { + x_high = x_new[i]; + x_low = 0.5 * (x_new[i-1] + x_new[i]); + } + else + { + x_high = 0.5 * (x_new[i] + x_new[i+1]); + x_low = 0.5 * (x_new[i-1] + x_new[i]); + } + + unsigned counter=0; + double sum = 0.0; + + + while(fft->get_spectrum(j, 0) >= x_low && fft->get_spectrum(j, 0) <= x_high && jget_spectrum(j, 1); + ++j; } - //std::cout << std::endl << " total: " << total_counter << std::endl; - //std::cout << "omega_max " << "old: " << fft->get_spectrum(N_old-1, 0) - // << " new: " << x_new[N_new-1] << std::endl; + + if (counter == 0) // this avoids a NaN if the bin does not contain any data + { + y_new[i] += 0.0; + } + else + { + y_new[i] += sum/counter; + } + + total_counter += counter; + + } } diff --git a/src/include/large_index_storage.hpp b/src/include/large_index_storage.hpp index 12e8468..f8408f9 100644 --- a/src/include/large_index_storage.hpp +++ b/src/include/large_index_storage.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,20 +19,15 @@ */ +#pragma once -#include -#include -//#include "detector.hpp" - -#ifndef LARGE_INDEX_STORAGE_RPAUSCH -#define LARGE_INDEX_STORAGE_RPAUSCH /** - * \brief storage class by Richard - * + * \brief storage class by Richard + * * this class provides a storage systems for large arrays for which - * there only a few non zero values, which are additionaly close together \n + * there only a few non zero values, which are additionally close together \n * usage: Large_index_storage n */ @@ -41,32 +36,28 @@ template class Large_index_storage { public: -// constructor, destructor: - //! \brief constructor - /*! @param N number of posible non zero values (just single segment) - @param start starting point of array f.e. x[start - x] x >0 --> is not - defined */ - inline Large_index_storage(const unsigned N, const unsigned start); - + // constructor, destructor: + //! \brief constructor + /*! @param N number of possible non zero values (just single segment) + @param start starting point of array f.e. x[start - x] x >0 --> is not + defined */ + inline Large_index_storage(const unsigned N, + const unsigned start); + inline ~Large_index_storage(); - //! \brief direct access to data, returns reference - /*! \param i is equivalent to x[i] */ + //! \brief direct access to data, returns reference + /*! \param i is equivalent to x[i] */ inline T& operator[](unsigned i); - //! \brief returns stored data and zero if not in [start, start + N -1] - /*! \param i is equivalent to x[i] */ + //! \brief returns stored data and zero if not in [start, start + N -1] + /*! \param i is equivalent to x[i] */ inline T operator()(unsigned i); private: - const unsigned N, start; - T* signal; - + const unsigned N, start; + T* signal; }; #include "large_index_storage.tpp" - - -#endif - diff --git a/src/include/large_index_storage.tpp b/src/include/large_index_storage.tpp index e11eecf..a16a263 100644 --- a/src/include/large_index_storage.tpp +++ b/src/include/large_index_storage.tpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -20,43 +20,41 @@ -#include "large_index_storage.hpp" - - - template< typename T> -inline Large_index_storage::Large_index_storage(const unsigned N, const unsigned start) - : N(N), start(start), signal(new T[N]) +inline Large_index_storage::Large_index_storage(const unsigned N, + const unsigned start) + : N(N), + start(start), + signal(new T[N]) { } - + template inline Large_index_storage:: ~Large_index_storage() { delete[] signal; } -//! \brief direct access to data, returns reference +//! \brief direct access to data, returns reference /*! \param i is equivalent to x[i] */ template inline T& Large_index_storage::operator[](unsigned i) { - if (!((i>=start) && (i " - << start << " - " << start+N<=start) && (i " + << start << " - " << start+N< inline T Large_index_storage::operator()(unsigned i) { if((i>= start) && (i< start+N)) return signal[i-start]; // returns data - else + else return T(0); // returns zero (or equivalent) } - - diff --git a/src/include/load_txt.hpp b/src/include/load_txt.hpp index af46135..fabe034 100644 --- a/src/include/load_txt.hpp +++ b/src/include/load_txt.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,75 +19,62 @@ */ +#pragma once -#include -#include -//#include - - -#ifndef LOAD_TXT_RPAUSCH -#define LOAD_TXT_RPAUSCH - -#include "import_from_file.hpp" +#include using namespace std; - //! \brief function returning the number of lines of a file /*! @param target a 'const char*' containing file location */ -unsigned linecounter(const char* target) +unsigned linecounter(const char* target) { - unsigned counter = 0; // linecounter - std::ifstream file(target); // target file - while (!file.eof()) { - if (file.get() == '\n') { ++counter; } - } - file.close(); // unnecessary but OK - return counter + 1; // add one for last line without '\n' + unsigned counter = 0; // line counter + std::ifstream file(target); // target file + while (!file.eof()) + { + if (file.get() == '\n') + ++counter; + } + file.close(); // unnecessary but OK + + return counter + 1; // add one for last line without '\n' } - - - -void load_txt( const char target[], const unsigned linenumber, one_line* data) -{ - ifstream file(target); // file to load - - string storage; // storage of short strings - for (unsigned i=0; !file.eof(); ++i) - { - // check for FORTRAN error: 1.234e-123 is stored wrongly as 1.234-123 ! -#ifndef CHECK_FOR_FORTRAN_ERROR - file >> data[i/7].intern_data[i%7]; +void load_txt( const char target[], + const unsigned linenumber, + one_line* data) +{ + ifstream file(target); // file to load + + string storage; // storage of short strings + for (unsigned i=0; !file.eof(); ++i) + { + // check for FORTRAN error: 1.234e-123 is stored wrongly as 1.234-123 ! +#ifndef CHECK_FOR_FORTRAN_ERROR + file >> data[i/7].intern_data[i%7]; #else - file >> storage; - // going through the string and checking if the Mathematica output - // error occurred: - for(unsigned j=1; storage[j] != '\0'; ++j){ // ignoring a first sign - if (storage[j] == '-' && storage[j-1] != 'e'){ - std::string str1, str2; - str1 = storage.substr(0, j); // string befor error - str2 = storage.substr(j); // string after error - storage = str1 + 'e' + str2; // corrected output - } - } - data[i/7].intern_data[i%7] = atof(storage.data()); - // storing data (7 doubles per line) (string to double) -#endif + file >> storage; + // going through the string and checking if the Mathematica output + // error occurred: + for(unsigned j=1; storage[j] != '\0'; ++j) + { // ignoring a first sign + if (storage[j] == '-' && storage[j-1] != 'e') + { + std::string str1, str2; + str1 = storage.substr(0, j); // string before error + str2 = storage.substr(j); // string after error + storage = str1 + 'e' + str2; // corrected output + } } - - - file.close(); - //std::cout << "Done\n\n"; - - ///////////////////////////////////////////////////////////// - // GOT DATA ///////////////////////////////////////////////// - ///////////////////////////////////////////////////////////// -} - - + data[i/7].intern_data[i%7] = atof(storage.data()); + // storing data (7 doubles per line) (string to double) #endif + } + file.close(); + +} diff --git a/src/include/makefile b/src/include/makefile index 53d3218..cfa61bc 100644 --- a/src/include/makefile +++ b/src/include/makefile @@ -1,4 +1,4 @@ -# Copyright 2014 Richard Pausch +# Copyright 2014-2016 Richard Pausch, Alexander Koehler # # This file is part of Clara 2. # @@ -17,20 +17,19 @@ # If not, see . # -# make the detectorclass correctly +# make the detector class correctly CC = g++ CFLAGS = -Wall -O3 -c +# compile against fftw library +# http://www.fftw.org/ - an open source FFT library under GPL 2.0 license +CFFT = -lfftw3 -lm +all: libDetector.a fileExists.o -all: libDetector.a - -libDetector.a: detector_e_field.o detector_dft.o detector_fft.o +libDetector.a: detector_e_field.o detector_dft.o detector_fft.o rm -f libDetector.a - ar cr libDetector.a detector_e_field.o detector_dft.o detector_fft.o - -#copy: libDetector.a -# cp libDetector.a .. + ar cr libDetector.a detector_e_field.o detector_dft.o detector_fft.o # E_field detector: @@ -42,10 +41,11 @@ detector_dft.o: detector_dft.cpp detector_dft.hpp vector.hpp utilities.hpp $(CC) $(CFLAGS) detector_dft.cpp # Fast Fourier Transformation detectors: -detector_fft.o: detector_fft.cpp detector_fft.hpp vector.hpp utilities.hpp fft_ned.hpp - $(CC) $(CFLAGS) detector_fft.cpp - +detector_fft.o: detector_fft.cpp detector_fft.hpp vector.hpp utilities.hpp ned_fft.hpp ../settings.hpp + $(CC) $(CFLAGS) $(CFFT) detector_fft.cpp +fileExists.o: fileExists.hpp fileExists.cpp + $(CC) $(CFLAGS) fileExists.cpp # How do I include the detector.hpp <-- vector.hpp dependency diff --git a/src/include/ned_fft.hpp b/src/include/ned_fft.hpp new file mode 100644 index 0000000..1df3cf5 --- /dev/null +++ b/src/include/ned_fft.hpp @@ -0,0 +1,218 @@ +/** + * Copyright 2014-2017 Richard Pausch + * + * This file is part of Clara 2. + * + * Clara 2 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 3 of the License, or + * (at your option) any later version. + * + * Clara 2 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 should have received a copy of the GNU General Public License + * along with Clara 2. + * If not, see . + */ + +#include +#pragma once + + + +inline unsigned power_of_two(unsigned N) +{ + unsigned exponent=1; + for(; N > (1u< // A...time, T...data +class ned_FFT +{ + +public: + // constructor + ned_FFT(unsigned N, + A x_original[], + T y_original[]) + : N_data(N), + x_equi(0), + y_equi(0), + data_complex(0), + spektrum(0), + omega(0) + { + x_equi = new A[N]; + y_equi = new T[N]; + + delta_t = interpolation_equi(x_original, y_original, N_data, + x_equi, y_equi, N_data); + + unsigned exponent=1; + for(; N_data > (1u< x_0[i]) + { + std::cerr << "error 01: interpolation inverted (ned_fft.hpp) " + << i << " --> " << x_0[i-1] + << " <=! " << x_0[i] << "\n"; + } + } + + const A min = x_0[0]; + const A max = x_0[N_0-1]; + + // creating equidistant x_values + for (unsigned i=0; i < N_1; ++i) + x_1[i] = min + (max-min)/(N_1) * i; + + // calculating y_values + unsigned j=0; + for (unsigned i=0; iA) + void spektrum_calc() + { + spektrum = new A[N_bin]; + for (unsigned i=0; i spektrum(N_data, x_data, y_data); diff --git a/src/include/physics_units.hpp b/src/include/physics_units.hpp index 35b96b2..df4cba5 100644 --- a/src/include/physics_units.hpp +++ b/src/include/physics_units.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -20,34 +20,27 @@ -// GLOBAL Variables +// GLOBAL Variables -#ifndef PHYSICS_UNITS_RPAUSCH -#define PHYSICS_UNITS_RPAUSCH +#pragma once namespace phy { -const double c = 299792458.; - /// speed of light [m/s] - -const double q = 1.602177E-19; - /// charge [ C = A*s ] + const double c = 299792458.; + /// speed of light [m/s] -/* -const double charge_e = 4.803204272E-10; - /// charge [ g^1/2 cm^3/2 s^-1 ] -*/ + const double q = 1.602177E-19; + /// charge [ C = A*s ] -const double m_e = 9.1093E-31; - /// mass of the electron [kg] - -const double epsilon_0 = 8.854188E-12; -/// + /* + const double charge_e = 4.803204272E-10; + /// charge [ g^1/2 cm^3/2 s^-1 ] + */ + const double m_e = 9.1093E-31; + /// mass of the electron [kg] + const double epsilon_0 = 8.854188E-12; + /// vacuum permittivity [F/m] } - - -#endif - diff --git a/src/include/string_manipulation.hpp b/src/include/string_manipulation.hpp deleted file mode 100644 index 74d71c8..0000000 --- a/src/include/string_manipulation.hpp +++ /dev/null @@ -1,64 +0,0 @@ -/** - * Copyright 2014 Richard Pausch - * - * This file is part of Clara 2. - * - * Clara 2 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 3 of the License, or - * (at your option) any later version. - * - * Clara 2 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 should have received a copy of the GNU General Public License - * along with Clara 2. - * If not, see . - */ - - - -#include -#include -#include -#include - -#ifndef formating_rpausch -#define formating_rpausch - -void formating(std::string& index_s, const std::string& digits_s) -{ - // ugly formating: - unsigned index_u = atoi(index_s.c_str()); // make a number - unsigned digits_u = atoi(digits_s.c_str());// make a number - - unsigned dummy=1; - for(unsigned i = 0; i< digits_u; ++i) - dummy *=10; - - if(index_u >= dummy) std::cout << " WARNING: index too large" << std::endl; - - std::string format = ("%0" + digits_s + "d"); // format style - char* index_c = new char[digits_u+5]; // create char dummy - sprintf(index_c, format.c_str(), index_u); // confert frum unsigned to char - index_s = index_c; // confert char to string - delete[] index_c; // delte char dummy - // end formating -} - - -double get_double(std::string txt) -{ - return atof(txt.c_str()); -} - - - -unsigned get_unsigned(std::string txt) -{ - return (unsigned)atoi(txt.c_str()); -} - -#endif diff --git a/src/include/utilities.hpp b/src/include/utilities.hpp index ecda4ff..14c41ed 100644 --- a/src/include/utilities.hpp +++ b/src/include/utilities.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,42 +19,35 @@ */ +#pragma once -#ifndef UTILITIES_RPAUSCH -#define UTILITIES_RPAUSCH - -namespace util { - - //goal: to increase readability of code - - template /// a generic square function - inline A square(A a ) - { - return a*a; - } - - template /// a more generic square function - inline R square(A a ) - { - return a*a; - } - - - template /// a generic cube function - inline A cube(A a) - { - return a*a*a; - } - - template /// a more generic cube function - inline R cube(A a) - { - return a*a*a; - } - -} +namespace util +{ + + //goal: to increase readability of code -#endif + template /// a generic square function + inline A square(A a ) + { + return a*a; + } + template /// a more generic square function + inline R square(A a ) + { + return a*a; + } + template /// a generic cube function + inline A cube(A a) + { + return a*a*a; + } + template /// a more generic cube function + inline R cube(A a) + { + return a*a*a; + } + +} diff --git a/src/include/vector.hpp b/src/include/vector.hpp index a5bc2b9..bbb35f4 100644 --- a/src/include/vector.hpp +++ b/src/include/vector.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,6 +19,7 @@ */ +#pragma once #include #include @@ -26,13 +27,11 @@ #include -#ifndef VECTOR_RPAUSCH -#define VECTOR_RPAUSCH /** - * \brief Vector class by Richard - * - * This class provied operations like: +, -, *, %, /, =, +=, [], mag() \n + * \brief Vector class by Richard + * + * This class provides operations like: +, -, *, %, /, =, +=, [], mag() \n * usage: Vector \n * short for Vector --> R_vec */ @@ -41,55 +40,57 @@ template< typename T, unsigned N> class Vector{ public: -// constructor, destructor: - /// basic constructor settting everything to zero - inline Vector(); - - //! \brief special constructor for 3D - /*! @param x first coordinate - @param y second coordinate - @param z third coordinate */ - inline Vector(double x, double y=0., double z=0.); - - /// copy constructor --> not neccesairy - - /// destructor - inline ~Vector() {} // not neccesary - - -// Getters and Setters: - /// access data - inline T & operator[](unsigned i); - - /// access data - const inline T & operator[](unsigned i) const; - - - -// Calculations: - - /// addition + // constructor, destructor: + /// basic constructor setting everything to zero + inline Vector(); + + //! \brief special constructor for 3D + /*! @param x first coordinate + @param y second coordinate + @param z third coordinate */ + inline Vector(double x, + double y=0., + double z=0.); + + /// copy constructor --> not necessary + + /// destructor + inline ~Vector() {} // not necessary + + + // Getters and Setters: + /// access data + inline T & operator[](unsigned i); + + /// access data + const inline T & operator[](unsigned i) const; + + + + // Calculations: + + /// addition inline Vector operator+ (const Vector& v) const; - /// subtraction + /// subtraction inline Vector operator- (const Vector& v) const; - - /// magnitude + + /// magnitude inline T mag() const; - /// scalar product + /// scalar product inline T dot(const Vector& v) const; - - /// multiplication with scalar + + /// multiplication with scalar inline Vector dot(const T a) const; - - /// cross-product-warning + + /// cross-product-warning inline Vector cross(const Vector& v) const; - - /// assign addition + + /// assign addition inline Vector & operator += (const Vector& v); - inline Vector unit_vec(); + inline Vector unit_vec() const; //////////// @@ -99,59 +100,63 @@ class Vector{ /// make complex vector real (absolute value of complex number) Vector abs() const; - - - + + + private: - /// stored data - T data[N]; + /// stored data + T data[N]; }; -// Member function spezialization: +// Member function specialization: /// 3D cross product template<> -inline Vector Vector::cross(const Vector& v) const; +inline Vector Vector::cross(const Vector& v) const; // global methods --> symbols for calculations and printing template< typename T, unsigned N> /// Vector * Vector --> scalar product -inline double operator * (const Vector& a , const Vector& b ); +inline double operator * (const Vector& a , + const Vector& b ); template< typename T, unsigned N> /// Vector * scalar -inline Vector operator * (const Vector & v, const double a); +inline Vector operator * (const Vector & v, + const double a); template< typename T, unsigned N> /// scalar * Vector -inline Vector operator * (const double a, const Vector & v); +inline Vector operator * (const double a, + const Vector & v); template< typename T, unsigned N> /// Vector / scalar -inline Vector operator / (const Vector & v, const double a); +inline Vector operator / (const Vector & v, + const double a); template< typename T, unsigned N> /// cross product --> Vector % Vector -inline Vector operator % (const Vector & a, const Vector & b); +inline Vector operator % (const Vector & a, + const Vector & b); //! \brief output stream used on vector object -/*! @param os output stream - @param v vector */ +/*! @param os output stream + @param v vector */ template< typename T, unsigned N> /// print Vector -inline std::ostream & operator << (std::ostream & os, const Vector & v); +inline std::ostream & operator << (std::ostream & os, + const Vector & v); /*! \var typedef Vector R_vec - \brief A type definition for a 3D double vector. - - Because in physics this is the most widly used vector, there is a special - typedef for it. - */ + \brief A type definition for a 3D double vector. + + Because in physics this is the most widely used vector, there is a special + typedef for it. +*/ #include "vector.tpp" typedef Vector R_vec; - -#endif diff --git a/src/include/vector.tpp b/src/include/vector.tpp index a302ed5..2f955c7 100644 --- a/src/include/vector.tpp +++ b/src/include/vector.tpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -27,16 +27,17 @@ //basic constructor template< typename T, unsigned N> -inline Vector::Vector() +inline Vector::Vector() { - for (unsigned i=0; i -inline Vector::Vector(double x, double y, double z) +inline Vector::Vector(double x, + double y, + double z) { assert( N==3 ); // better via template? one more dummy function data[0]=x; @@ -44,7 +45,7 @@ inline Vector::Vector(double x, double y, double z) data[2]=z; } - + //////////////////////////////////////////////// @@ -52,15 +53,15 @@ inline Vector::Vector(double x, double y, double z) // access data template< typename T, unsigned N> -inline T & Vector::operator[](unsigned i) +inline T & Vector::operator[](unsigned i) { assert(i -const inline T & Vector::operator[](unsigned i) const +template< typename T, unsigned N> +const inline T & Vector::operator[](unsigned i) const { assert(i::operator[](unsigned i) const -///////////////////////////////////////////////////// +///////////////////////////////////////////////////// // Calculations: - + /// addition -template< typename T, unsigned N> -inline Vector Vector::operator+ (const Vector& v) const +template< typename T, unsigned N> +inline Vector Vector::operator+ (const Vector& v) const { Vector back; - for (unsigned i=0; i -inline Vector Vector::operator- (const Vector& v) const +template< typename T, unsigned N> +inline Vector Vector::operator- (const Vector& v) const { Vector back; - for (unsigned i=0; i -inline T Vector::mag() const -{ +template< typename T, unsigned N> +inline T Vector::mag() const +{ T sqr_magnitude=0; - for (unsigned i=0; i -inline T Vector::dot(const Vector& v) const +template< typename T, unsigned N> +inline T Vector::dot(const Vector& v) const { T result=0; - for (unsigned i=0; i -inline Vector Vector::dot(const T a) const +template< typename T, unsigned N> +inline Vector Vector::dot(const T a) const { Vector result; - for (unsigned i=0; i -inline Vector Vector::cross(const Vector& v) const +template< typename T, unsigned N> +inline Vector Vector::cross(const Vector& v) const { std::cout << std::endl << "---------------------------" << std::endl << @@ -143,22 +144,22 @@ inline Vector Vector::cross(const Vector& v) const Vector dummy; return dummy; } - + /// assign addition -template< typename T, unsigned N> -inline Vector & Vector::operator += (const Vector& v) +template< typename T, unsigned N> +inline Vector & Vector::operator += (const Vector& v) { - for (unsigned i=0; i -inline Vector Vector::unit_vec() -{ +template< typename T, unsigned N> +inline Vector Vector::unit_vec() const +{ return *this / mag() ; } @@ -166,7 +167,7 @@ inline Vector Vector::unit_vec() // make real vector complex -template< typename T, unsigned N> +template< typename T, unsigned N> Vector, N> Vector::make_complex() const { Vector, N> result; @@ -177,27 +178,27 @@ Vector, N> Vector::make_complex() const } // make complex vector real -template< typename T, unsigned N> +template< typename T, unsigned N> Vector Vector::abs() const { Vector result; for(unsigned i = 0; i< N; ++i) result[i] = std::abs(data[i]); - + return result; } /// 3D cross product template<> -inline Vector Vector::cross(const Vector& v) const +inline Vector Vector::cross(const Vector& v) const { - Vector result; - result[0] = data[1]*v[2]-v[1]*data[2]; - result[1] = data[2]*v[0]-v[2]*data[0]; - result[2] = data[0]*v[1]-v[0]*data[1]; - - return result; + Vector result; + result[0] = data[1]*v[2]-v[1]*data[2]; + result[1] = data[2]*v[0]-v[2]*data[0]; + result[2] = data[0]*v[1]-v[0]*data[1]; + + return result; } @@ -206,47 +207,53 @@ inline Vector Vector::cross(const Vector& v) co template< typename T, unsigned N> /// Vector * Vector -inline double operator * (const Vector& a , const Vector& b ) /// Skalarprodukt +inline double operator * (const Vector& a , + const Vector& b ) /// scalar product { - return a.dot(b); + return a.dot(b); } template< typename T, unsigned N> /// Vector * scalar -inline Vector operator * (const Vector & v, const double a) +inline Vector operator * (const Vector & v, + const double a) { - return v.dot(a); + return v.dot(a); } template< typename T, unsigned N> /// scalar * Vector -inline Vector operator * (const double a, const Vector & v) +inline Vector operator * (const double a, + const Vector & v) { - return v.dot(a); + return v.dot(a); } template< typename T, unsigned N> /// Vector / scalar -inline Vector operator / (const Vector & v, const double a) +inline Vector operator / (const Vector & v, + const double a) { - return v.dot(1/a); + return v.dot(1/a); } template< typename T, unsigned N> /// cross product --> Vector % Vector -inline Vector operator % (const Vector & a, const Vector & b) // cross-product +inline Vector operator % (const Vector & a, + const Vector & b) // cross-product { - return a.cross(b); + return a.cross(b); } //! \brief output stream used on vector object -/*! @param os output stream - @param v vector */ +/*! @param os output stream + @param v vector */ template< typename T, unsigned N> /// print Vector -inline std::ostream & operator << (std::ostream & os, const Vector & v) +inline std::ostream & operator << (std::ostream & os, + const Vector & v) { - os << "(" ; - for (unsigned i=0; i<(N-1); ++i) { - os << v[i] << " , "; - } - os << v[N-1] << ")"; - return os; + os << "(" ; + for (unsigned i=0; i<(N-1); ++i) + { + os << v[i] << " , "; + } + os << v[N-1] << ")"; + return os; } - diff --git a/src/main.cpp b/src/main.cpp index a661ace..ca83a85 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,14 +19,13 @@ */ - - #include #include #include "parallel_jobs.h" #include "all_directions.hpp" + int main(void) { printf("start\n"); @@ -35,7 +34,7 @@ int main(void) /** the function 'start_array()' calls the parallelization procedure * defined in 'parallel_jobs.h'. The variables 'numtasks' and 'rank' - * are set by this function and allow to use the parallel distrubuted + * are set by this function and allow to use the parallel distributed * indexes to be used. **/ if(start_array(&numtasks, &rank) != 0) @@ -44,9 +43,9 @@ int main(void) ////////////////////////////////////////////////// // do some work here ////////////////////////////////////////////////// - + /* hostname of the machine this task is running on */ - char * pHost; + char * pHost; pHost = getenv("MYHOSTNAME"); @@ -61,60 +60,62 @@ int main(void) const unsigned int N_max = 2001; /* max number of trajectories */ unsigned int i; /* index */ for(i=rank; i MPI + * __PARALLEL_SETTING__ == 2 --> PBS array job + */ #if __PARALLEL_SETTING__ == 1 +/* only include mpi header only if MPI was selected */ #include "mpi.h" #endif -int start_array(int* numtasks, int* rank) +/** this function emulates a PBS Array job task and either runs it + * in a PBS array environment or in an MPI environment + * + * @param numtasks pointer to int where number of total parallel + * tasks should be stored + * @param rank pointer to int where the rank/task-id of a specific + * task should be stored + * @return int with error code: 0 - >successful + * 1 -> MPI error + * 2 -> numtask not allocated + * 3 -> rank not allocated + */ +int start_array(int* numtasks, + int* rank) { + /* check if numtask is a valid pointer */ if(!numtasks) return 2; + /* check if rank is a valid pointer */ if(!rank) return 3; +/* in case MPI is used for parallelization */ #if __PARALLEL_SETTING__ == 1 - // MPI used for parallization + /* init MPI on this core */ int rc = MPI_Init(NULL, NULL); + /* check if MPI init was successful */ if (rc != MPI_SUCCESS) - { - printf("Error starting MPI program. Terminating program!\n"); - MPI_Abort(MPI_COMM_WORLD, rc); - return 1; - } + { + printf("Error starting MPI program. Terminating program!\n"); + MPI_Abort(MPI_COMM_WORLD, rc); + return 1; + } + /* get MPI numtasks and rank */ MPI_Comm_size(MPI_COMM_WORLD, numtasks); MPI_Comm_rank(MPI_COMM_WORLD, rank); +/*in case PBS array jobs are used for parallelization */ #elif __PARALLEL_SETTING__ == 2 - // Array jobs used for parallization - char* dump; + char* dump; /* temporary memory pointer for reading bash + environment variables */ + /* get rank = PBS_ARRAYID */ dump = getenv("PBS_ARRAYID"); *rank = atoi(dump); + /* get numtask = ARRAYMAX */ dump = getenv("ARRAYMAX"); *numtasks = atoi(dump); +/* throw compile time error if neither MPI nor PBS array jobs are selected */ #else - // non of the above selected #error parallel setting not suported #endif @@ -67,38 +100,42 @@ int start_array(int* numtasks, int* rank) return 0; } + +/** function to clean up after parallel job (if needed) + * @return int error code: 0 -> successful + */ int end_array(void) { +/* in case MPI is used for parallelization */ #if __PARALLEL_SETTING__ == 1 - // MPI used for parallization - MPI_Finalize(); - + MPI_Finalize(); +/* in case PBS array jobs are used for paralellization */ #elif __PARALLEL_SETTING__ == 2 - // Array jobs used for parallization - + /* nothing needs to be done */ +/* throw compile time error if neither MPI nor PBS array jobs are selected */ #else - // non of the above selected #error parallel setting not suported #endif - return 0; + return 0; } -/* + +/** function to check if a specific file is found in the working + * directory that forces a soft stop + * + * @return int 0 -> file not found (continue) + * 1 -> file found (stop running) + */ int check_break(void) { - char stop[] = "break.now"; - FILE* file = fopen(stop, "r"); - if(file != 0) - { - fclose(file); - //perror("breakpoint found \n"); - return 1; - } - - return 0; + char stop[] = "break.now"; /* file name that causes a stop */ + FILE* file = fopen(stop, "r"); /* try read access to file */ + if(file != 0) /* file found */ + { + fclose(file); /* close file handler again */ + return 1; /* return 1 = file found */ + } + + return 0; /* file not found - return 0 */ } -*/ - -#endif - diff --git a/src/prepare_job.sh b/src/prepare_job.sh index eec6520..adc665c 100755 --- a/src/prepare_job.sh +++ b/src/prepare_job.sh @@ -1,6 +1,6 @@ #!/bin/bash -# Copyright 2014 Richard Pausch +# Copyright 2014-2016 Richard Pausch # # This file is part of Clara 2. # @@ -24,26 +24,47 @@ #--------------------------------------------------------------- # INFO: # Running this shell script generates a submit file either -# for a MPI parallelsation or a PBS-Array-Job parallelsation. +# for a MPI parallelization or a PBS-Array-Job parallelization. #--------------------------------------------------------------- +#------------ +# USER SETUP: +#------------ +# name of cluster job: JOBNAME=MPI_lib1 + +# name of submit file: SUBMITFILE=submit.sh + +# name of queue to submit job to: QUEUE=laser + +# name of file with environment setup: MODULES2LOAD=clara2_hypnos.modules +# maximum simulation time WALLTIME=02:00:00 -# Jobstructure: -# for MPI: +# Job parallel structure: +# for MPI (number of node with (each) number of cores): NUMBERNODES=4 NUMBERCORES=32 # for Array Jobs: +# set number of parallel tasks to the same as with MPI ARRAYMAX=$(echo $NUMBERNODES " * " $NUMBERCORES " - 1 " | bc ) -# User interaction required +#--------------- +# END USER SETUP +#--------------- + + +#------------------------- +# RUN TIME USE INTERACTION +#------------------------- + +# user interaction required to select parallelization # present options: echo "for MPI enter 1" echo "for PBS-Array jobs enter 2" @@ -57,11 +78,16 @@ case $DECISION in exit 1 ;; esac +#-------------------- +# END USE INTERACTION +#-------------------- -source $MODULES2LOAD +#----------------------- +# GENERATE SUBMIT SCRIPT +#----------------------- -# generate submit file (the same for MPI and PBS-array jobs) +# generate submit file (here the same for MPI and PBS-array jobs) echo "#PBS -q "$QUEUE > $SUBMITFILE echo "#PBS -l walltime="$WALLTIME >> $SUBMITFILE echo "#PBS -N "$JOBNAME >> $SUBMITFILE @@ -71,7 +97,7 @@ echo "#PBS -e ./error.txt" >> $SUBMITFILE echo "#PBS -o ./output.txt" >> $SUBMITFILE echo "" >> $SUBMITFILE -# in case of a MPI parallelisation: +# in case of a MPI parallelization: if [ $DECISION -eq 1 ] then echo "#PBS -l nodes="$NUMBERNODES":ppn="$NUMBERCORES >> $SUBMITFILE @@ -93,18 +119,35 @@ if [ $DECISION -eq 1 ] then echo -n "mpiexec --prefix $MPIHOME -tag-output --display-map -x LIBRARY_PATH -x LD_LIBRARY_PATH -n " >> $SUBMITFILE echo $(echo $NUMBERNODES " * " $NUMBERCORES | bc )" ./executable | tee output.log" >> $SUBMITFILE - make - make MPI else - # in case of PBS-array paparllelisation - echo " ./executable" >> $SUBMITFILE - make - make ARRAY + # in case of PBS-array parallelization + echo " ./executable" >> $SUBMITFILE fi # add newline to submit file echo "" >> $SUBMITFILE +#--------------------------- +# END GENERATE SUBMIT SCRIPT +#--------------------------- +#------------- +# COMPILE CODE +#------------- + +# load bash environment for compilation +source $MODULES2LOAD + +# ISSUE #65 why is make required twice? +make +if [ $DECISION -eq 1 ] +then + make MPI +else + make ARRAY +fi +#----------------- +# END COMPILE CODE +#----------------- diff --git a/src/process_data.cpp b/src/process_data.cpp index 29b9732..3065f22 100644 --- a/src/process_data.cpp +++ b/src/process_data.cpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch, Alexander Koehler * * This file is part of Clara 2. * @@ -19,194 +19,203 @@ */ - -#include -#include #include -#include -#include "gzip_lib.hpp" +#include "include/input_output.hpp" +#include "settings.hpp" +#include "setFilename.hpp" + +/** data structure to hold spectral data + * that is automatically initialized with zeros + */ template struct spectrum { + /** constructor + * initialize data with zeros + */ spectrum(void) { + /* initialize all spectral data with zeros */ for(unsigned int i=0; i successful + */ +int main() { - - if(argc != 2) - { - std::cout << "wrong usage: needs 1 parameters not " << argc -1 << std::endl; - std::cout << "give: encoding" << std::endl; - return 1; - } - - using namespace std; - const unsigned int index_files_first = 0; - const unsigned int index_files_last = 2000; - const unsigned int N_omega = 2048; - const unsigned int N_theta = 120; - const unsigned int N_phi = 2; - const unsigned int N_direction = N_theta*N_phi; - const char input_pattern[] = "my_spectrum_trace%06d.dat"; - const char output_pattern[] = "my_spectrum_all_%03d.dat"; - const unsigned int N_split = 8; - - spectrum* data = new spectrum[N_direction]; + /* compute total number of observation direction */ + const unsigned int N_direction = param::N_theta * param::N_phi; + /* spectral container to hold total spectrum (for each direction) */ + spectrum* data = new spectrum[N_direction]; + /* ----------------------------- + READ SPECTRA FROM INPUT FILES + -----------------------------*/ + for(unsigned int index_files = param::index_files_first; index_files <= param::index_files_last; ++index_files) + { + /* set file name for spectra from trace with id index_files */ + char filename[param::N_char_filename]; + setFilename(filename, param::outputFileTemplate, index_files, param::N_char_filename); - // -------- get store info ----------- - bool asci_input; - std::string store_str = argv[1]; - if(!store_str.compare("asci")) + /* ----- read input file ------ */ + if(param::ascii_output) /* if files loaded have ASCII format */ { - asci_input = true; - std::cout << "ASCI input" << std::endl; + FILE* pFile = fopen(filename, "r"); /* file handler */ + + if(pFile == NULL) /* file handler doe not point to file - abort this file */ + { + std::cout << "abort ascii file " << index_files + << "/" << param::index_files_last << std::endl; + continue; + } + else /* file handler points to valid file */ + { + std::cout << "load ascii file " << index_files + << "/" << param::index_files_last << std::endl; + } + + /* process all directions */ + for(unsigned int j=0; j phi-index:" << index_phi << std::endl; - throw "error output"; - } + std::cerr << "error: could not write output-file --> phi-index:" << index_phi << std::endl; + throw "error output"; } - - - + } delete[] data; + /* -------------------------------------- + END STORE INCOHERENT RADIATION SPECTRA + --------------------------------------*/ - std::cout << "deleting files" << std::endl; - for(unsigned int index_files = index_files_first; index_files <= index_files_last; ++index_files) - { - char filename[256]; - sprintf(filename, input_pattern, index_files); + /* ------------------------------------- + REMOVE SPECTRA FROM INDIVIDUAL TRACES + -------------------------------------*/ - if(remove(filename) == 0) - std::cout << "removed: " << filename << std::endl; - else - std::cerr << "error removing: " << filename << std::endl; - } + std::cout << "deleting files" << std::endl; /* verbose output */ + /* for each file previously loaded: */ + for(unsigned int index_files = param::index_files_first; + index_files <= param::index_files_last; + ++index_files) + { + /* set file name */ + char filename[param::N_char_filename]; + setFilename(filename, param::outputFileTemplate, index_files, param::N_char_filename); + + /* try deleting and give verbose output about success or failure */ + if(remove(filename) == 0) + std::cout << "removed: " << filename << std::endl; + else + std::cerr << "error removing: " << filename << std::endl; + } + /* ----------------------------------------- + END REMOVE SPECTRA FROM INDIVIDUAL TRACES + -----------------------------------------*/ return 0; } diff --git a/src/run_through_data.hpp b/src/run_through_data.hpp index 6a9a3eb..f30441d 100644 --- a/src/run_through_data.hpp +++ b/src/run_through_data.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,127 +19,86 @@ */ - - -#include -#include - - -#ifndef RUN_THROUGH_DATA_RPAUSCH -#define RUN_THROUGH_DATA_RPAUSCH - +#pragma once #include "discrete.hpp" -#include "import_from_file.hpp" -#include "physics_units.hpp" #include "vector.hpp" +#include "settings.hpp" + +/** function calculating beta * gamma from beta + * + * @param R_vec beta (speed normalized to the speed of light + * @retun R_vec beta times relativistic gamma factor + * = momentum / (mass * speed of light) + */ inline R_vec beta_times_gamma(R_vec beta) { return std::sqrt(1./(1.-util::square(beta))) * beta; } - - - - +/** function to step through loaded data and add radiation amplitude to + * a given detector + * + * @param data (one_line type pointer) pointing to trajectory data + * @param linenumber number of time steps in data + * @param detector Clara detector for one specific direction + */ template -void run_through_data(const one_line* data, const unsigned linenumber, const unsigned N_angle, - DET detector) +void run_through_data(const one_line* data, + const unsigned linenumber, + DET detector) { - /* ---------- storing data : comparable to real data stream (not loke a file here) --- */ - - // USING: SI-UNITS !!! - - //Zeiten s --> s - Discrete time_fill(data[0].intern_data[6] /* *1.E-15 */ , - data[1].intern_data[6] /* *1.E-15 */ , - data[2].intern_data[6] /* *1.E-15 */ , - data[3].intern_data[6] /* *1.E-15 */ ); - - //Ort: in m --> m - Discrete location( R_vec(data[0].intern_data[0] /* *1.E-2 */ , - data[0].intern_data[1] /* *1.E-2 */ , - data[0].intern_data[2] /* *1.E-2 */ ), - R_vec(data[1].intern_data[0] /* *1.E-2 */ , - data[1].intern_data[1] /* *1.E-2 */ , - data[1].intern_data[2] /* *1.E-2 */ ), - R_vec(data[2].intern_data[0] /* *1.E-2 */ , - data[2].intern_data[1] /* *1.E-2 */ , - data[2].intern_data[2] /* *1.E-2 */ ), - R_vec(data[3].intern_data[0] /* *1.E-2 */ , - data[3].intern_data[1] /* *1.E-2 */ , - data[3].intern_data[2] /* *1.E-2 */ ), - &time_fill); - - //Impuls: beta*gamma --> phy::m_e*beta_times_gamma(beta)*phy:c --> kg*m/s - double btom = phy::m_e*phy::c; - Discrete momentum( btom*beta_times_gamma(R_vec(data[0].intern_data[3], - data[0].intern_data[4], - data[0].intern_data[5]) ), - btom*beta_times_gamma(R_vec(data[1].intern_data[3], - data[1].intern_data[4], - data[1].intern_data[5]) ), - btom*beta_times_gamma(R_vec(data[2].intern_data[3], - data[2].intern_data[4], - data[2].intern_data[5]) ), - btom*beta_times_gamma(R_vec(data[3].intern_data[3], - data[3].intern_data[4], - data[3].intern_data[5]) ), - &time_fill); - - More_discrete auto_fill(&time_fill); - - - Discrete beta(&time_fill); - Discrete gamma(&time_fill); - - gamma = auto_fill.momentum_to_gamma(momentum); - beta = auto_fill.momentum_to_beta(momentum, gamma); - - - - - /* -------- streaming the data and sending it to the detectors: ---------- */ - - - //std::cout << "running: "; - for(unsigned i=4; i time_fill; + const Discrete *time_fill_ref = &time_fill; + /* position */ + Discrete location( time_fill_ref ); + /* momentum */ + Discrete momentum( time_fill_ref ); + + /* decrived quantities */ + More_discrete auto_fill( time_fill_ref ); + /* beta */ + Discrete beta( time_fill_ref ); + /* gamma */ + Discrete gamma( time_fill_ref ); + + + namespace in = param::input; + /* step through data for each time step: streaming approach used (future option?) */ + for(unsigned i=0; i phy::m_e*beta_times_gamma(beta)*phy:c --> kg*m/s + const double btom = phy::m_e*phy::c; + momentum.next( btom*beta_times_gamma(R_vec(data[i].intern_data[in::index_beta_x] * in::convert_beta , + data[i].intern_data[in::index_beta_y] * in::convert_beta , + data[i].intern_data[in::index_beta_z] * in::convert_beta ) )); + + gamma.next(auto_fill.gamma(momentum.get_future())); + beta.next(auto_fill.beta(momentum.get_future(), gamma.get_future())); + + /* if all 4 entries of the Discrete class are filed with values + the radiation amplitude can be calculated */ + if(i >=3 ) { - for(unsigned k=0; kadd_to_spectrum(location.get_old(), + beta.get_old(), + beta.dot_old(), + time_fill.get_old(), + time_fill.get_delta_old()); } - //std::cout << std::endl; - -} - - -#endif + } +} diff --git a/src/setFilename.hpp b/src/setFilename.hpp new file mode 100644 index 0000000..d0fbf35 --- /dev/null +++ b/src/setFilename.hpp @@ -0,0 +1,48 @@ +/** + * Copyright 2016 Richard Pausch + * + * This file is part of Clara 2. + * + * Clara 2 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 3 of the License, or + * (at your option) any later version. + * + * Clara 2 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 should have received a copy of the GNU General Public License + * along with Clara 2. + * If not, see . + */ + +#pragma once + +#include "settings.hpp" + + +/** function to set file name with index based on template + * + * @param filename char pointer to char array where file name + * should be stored + * @param templateString pointer to char array with template + * @param index unsigned int with id number + * @param N_char number of characters available in filename array + */ +void setFilename(char* filename, + const char* templateString, + const unsigned int index, + const unsigned int N_char_filename) +{ + if(sprintf(filename, + templateString, + index) + > int(N_char_filename)) /* check if name fits in filename array */ + { + /* throw warning when buffer is to small for path name */ + std::cerr << "filename buffer too small!!! " << std::endl; + throw "Buffer to small!"; + } +} diff --git a/src/settings.hpp b/src/settings.hpp new file mode 100644 index 0000000..472703d --- /dev/null +++ b/src/settings.hpp @@ -0,0 +1,68 @@ +/** + * Copyright 2016 Alexander Koehler, Richard Pausch + * + * This file is part of Clara 2. + * + * Clara 2 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 3 of the License, or + * (at your option) any later version. + * + * Clara 2 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 should have received a copy of the GNU General Public License + * along with Clara 2. + * If not, see . + */ + + +#pragma once + +namespace param +{ +const double omega_max = 3.0e19; /* maximum of plotted frequency Hz */ +const double theta_max = 1.14594939; /* maximum of theta in degree */ +const unsigned int N_spectrum = 2048; /* number of frequencies "omega" */ +const unsigned int N_theta = 120; /* number of directions in first angle "theta" */ +const unsigned int N_phi = 2; /* number of directions in second angle "phi" */ +const unsigned int N_trace = 2000; /* maximum number of traces */ + +const unsigned int fft_length_factor = 1; /* needs to be a power of two */ + +const bool ascii_output = false; /* output spectra as ascii text */ + +const unsigned int N_char_filename=256; // number of characters +/* template for input traces (using C formatting): */ +const char traceFileTemplate[] = "/net/cns/projects/HPLsim/xray/debus/ELBEThomson/basicRun2/trace_%04d.txt"; +/* template for output spectra for each trace (using C formatting): */ +const char outputFileTemplate[] = "my_spectrum_trace%06d.dat"; + +namespace input +{ + const unsigned int index_time = 6; /* column id of time */ + + const unsigned int index_pos_x = 0; /* column id of x */ + const unsigned int index_pos_y = 1; /* column id of y */ + const unsigned int index_pos_z = 2; /* column id of z */ + + const unsigned int index_beta_x = 3; /* column id of beta_x */ + const unsigned int index_beta_y = 4; /* column id of beta_y */ + const unsigned int index_beta_z = 5; /* column id of beta_z */ + + /* unit conversion to SI units */ + const double convert_time = 1.0; /* s -> s */ + const double convert_pos = 1.0; /* m -> m */ + const double convert_beta = 1.0; /* None -> None */ +} + +// explicit for process_data +const unsigned int N_omega = N_spectrum; /* number of frequencies */ +const unsigned int index_files_first = 0; /* start index for reading trajectory spectra */ +const unsigned int index_files_last = N_trace; /* final index for reading trajectory spectra */ +/* template for output total (incoherent) spectra (using C formatting for phi index): */ +const char output_pattern[] = "my_spectrum_all_%03d.dat"; + +} // end namespace param diff --git a/src/single_direction.cpp b/src/single_direction.cpp new file mode 100644 index 0000000..332284b --- /dev/null +++ b/src/single_direction.cpp @@ -0,0 +1,144 @@ +/** + * Copyright 2014-2016 Richard Pausch + * + * This file is part of Clara 2. + * + * Clara 2 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 3 of the License, or + * (at your option) any later version. + * + * Clara 2 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 should have received a copy of the GNU General Public License + * along with Clara 2. + * If not, see . + */ + + +#include "single_direction.hpp" + + +#include "vector.hpp" +//#include "detector_e_field.hpp" /* currently not used */ +#include "detector_dft.hpp" +#include "detector_fft.hpp" + +/* only needed for near field calculation + * REMOVE ? */ +/* #include"large_index_storage.hpp" */ + +#include "interpolation.hpp" +#include "run_through_data.hpp" + + +/** + * calculates a single spectra for only one trace and one direction + * + * @param data pointer to trajectory data + * @param linenumber number of data points + * @param all_omega pointer to frequency values + * @param all_spectrum pointer to memory for spectra + * @param N_all_spec maximum number of spectral data allocated + * @param theta_offset offset of angle theta (used to set angle) + * @param phi_offset offset of angle phi (used to set angle) + **/ +int single_direction(const one_line* data, + const unsigned int linenumber, + const double* all_omega, + double* all_spectrum, + const unsigned N_all_spec, + const double theta_offset, + const double phi_offset) +{ + + /* ----------------------- detectors ------------------------ */ + /* creating the "looking vector" for the "detector class object" + * that "observes" the emitted radiation in the far field */ + + /* setup observation direction: */ + R_vec looking_vector; + double angle_theta; + double angle_phi; + + angle_theta = theta_offset/180.0 * M_PI; + + angle_phi = phi_offset/180.0 * M_PI; + + looking_vector = R_vec(std::cos(angle_theta) , + std::sin(angle_theta) * std::cos(angle_phi) , + std::sin(angle_theta) * std::sin(angle_phi) ); + + + /* -------- FFT ------------- */ + /* this performs the radiation calculation + * (for a single trace and a single direction) + * using the FFT algorithm */ + + + /* create memory for detectors */ + Detector_fft* detector_fft; + + /* create FFT detectors */ + detector_fft = new Detector_fft(looking_vector, linenumber); + + + + /* convert trajectory data to meaning full values */ + /* TO DO: FUNCTION CALLED FOR EACH DIRECTION - ISSUE #14 */ + run_through_data(data, linenumber, detector_fft); + + detector_fft->calc_spectrum(); + + + + /* this switches off the computation using a DFT method - ISSUE #13 */ +#if 0 + + /* ------------------ DFT ------------------- */ + /* right now, this is just a debugging method for the FFT + * it calculates the radiation for the exact same frequencies + * as the FFT */ + + + /* create detector objects for radiation */ + Detector_dft* detector_dft; + + + /* create DFT detectors (using omega[i] from FFT) to be on the exact same frequency */ + /* number of (use-full) frequencies from FFT */ + unsigned dummy_N = detector_fft->half_frequency(); + + /* get frequencies from FFT */ + double* dummy_omega = detector_fft->frequency; + detector_dft = new Detector_dft(looking_vector, dummy_N, dummy_omega); + + /* convert trajectory data to meaningful values */ + /* TO DO: FUNCTION CALLED FOR EACH DIRECTION - ISSUE #14 */ + run_through_data(data, linenumber, detector_dft); + + /* calculate spectra using DFT method */ + detector_dft->calc_spectrum(); + + +#endif + /* end DFT calculation switch - ISSUE #13 */ + + + + /* ------ interpolate spectra onto requested frequencies ---------- */ + /* this is the "integrated interpolation" found in interpolate.tpp */ + interpolation_int(detector_fft, all_omega, all_spectrum, N_all_spec); + + + /* ----- free memory used for detectors ---- */ + + /* TO DO: DFT detectors not used - ISSUE #13 */ + /* delete detector_dft; */ + delete detector_fft; + + return 0; +} diff --git a/src/single_trace.hpp b/src/single_direction.hpp similarity index 51% rename from src/single_trace.hpp rename to src/single_direction.hpp index e28df3f..c41b2c9 100644 --- a/src/single_trace.hpp +++ b/src/single_direction.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Richard Pausch + * Copyright 2014-2016 Richard Pausch * * This file is part of Clara 2. * @@ -19,39 +19,15 @@ */ - - -#include -#include -#include -#include - - - -#ifndef SINGLE_TRACE_RPAUSCH -#define SINGLE_TRACE_RPAUSCH - -using namespace std; - -#include "vector.hpp" - -#include "detector_e_field.hpp" -#include "detector_dft.hpp" -#include "detector_fft.hpp" +#pragma once #include "import_from_file.hpp" -/* only needed for near field calculation - * REMOVE ? */ -/* #include"large_index_storage.hpp" */ - -#include "physics_units.hpp" - /** * calculates a single spectra for only one trace and one direction * - * @param one_line pointer to trajectory data + * @param data pointer to trajectory data * @param linenumber number of data points * @param all_omega pointer to frequency values * @param all_spectrum pointer to memory for spectra @@ -59,23 +35,11 @@ using namespace std; * @param theta_offset offset of angle theta (used to set angle) * @param phi_offset offset of angle phi (used to set angle) **/ -int single_trace(const one_line* data, - const unsigned int linenumber, - const double* all_omega, - double* all_spectrum, - const unsigned N_all_spec, - const double theta_offset = 0.0, - const double phi_offset = 0.0); - - - -/* TO DO: this should be in a separate file - ISSUE #15 */ -/** - * check whether a file exists or not - * - * @param filename pointer to array containing file location - * @return Returs true if file exists, otherwise false. - **/ -bool file_exists(const char *filename); -#endif +int single_direction(const one_line* data, + const unsigned int linenumber, + const double* all_omega, + double* all_spectrum, + const unsigned N_all_spec, + const double theta_offset = 0.0, + const double phi_offset = 0.0); diff --git a/src/single_trace.cpp b/src/single_trace.cpp deleted file mode 100644 index 7d06c71..0000000 --- a/src/single_trace.cpp +++ /dev/null @@ -1,213 +0,0 @@ -/** - * Copyright 2014 Richard Pausch - * - * This file is part of Clara 2. - * - * Clara 2 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 3 of the License, or - * (at your option) any later version. - * - * Clara 2 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 should have received a copy of the GNU General Public License - * along with Clara 2. - * If not, see . - */ - - - - -#include "single_trace.hpp" - -#include "interpolation.hpp" -#include "discrete.hpp" -#include "run_through_data.hpp" - - - -/** - * calculates a single spectra for only one trace and one direction - * - * @param one_line pointer to trajectory data - * @param linenumber number of data points - * @param all_omega pointer to frequency values - * @param all_spectrum pointer to memory for spectra - * @param N_all_spec maximum number of spectral data allocated - * @param theta_offset offset of angle theta (used to set angle) - * @param phi_offset offset of angle phi (used to set angle) - **/ -int single_trace(const one_line* data, - const unsigned int linenumber, - const double* all_omega, - double* all_spectrum, - const unsigned N_all_spec, - const double theta_offset, - const double phi_offset) -{ - //////////////////////////////////////////////////////////// - /* CHANGE THIS: */ - /* this are the remainings of the old program - * structure, where several directions were calculated - * right here - ISSUE #12 */ - - const unsigned N_angle_theta = 1; - const unsigned N_angle_phi = 1; - - - ////////////////////////////////////////////////////////////// - - - - - /* ---------- time steps : delta_t ---------------------- */ - - /* get time step width "delta_t=stepwidth" from data */ - /* WARNING: HERE KNOWLEDGE ON DATA STRUCTURE IS REQUIRED - ISSUE #11 */ - double stepwidth = (data[6].intern_data[6] -data[5].intern_data[6]); - - - - - /* ----------------------- detectors ------------------------ */ - /* creating the "looking vector" for the "detector class object" - * that "observes" the emitted radiation in the far field */ - - /* creating different angles: */ - /* TO DO: SINCE ONLY A SINGLE DIRECION IS CALCULATED IN THIS FUNCTION - * THE ARRAY SETUP IS OBSOLETE - ISSUE #12 */ - R_vec looking_vector[N_angle_theta * N_angle_phi]; - double angle_theta[N_angle_theta]; - double angle_phi[N_angle_phi]; - - for(unsigned i=0; i< N_angle_theta; ++i) - { - angle_theta[i] = (double(i)/double(N_angle_theta)) * M_PI * 0.5 + theta_offset/180.0 * M_PI; - } - - for(unsigned i=0; i< N_angle_phi; ++i) - { - angle_phi[i] = (double(i)/double(N_angle_phi))*1.0 * M_PI + phi_offset/180.0 * M_PI; - } - - for (unsigned a=0; a