Hi,
Brief question: I have code that compiles without errors or warnings (using g++ on a Fedora machine in work), however when i run it there is segmentation fault at the following line:
cout << "Please input the minimum wavelength, in microns, for use in the Sellmeier Equation: " << endl;
I have moved the end of the program (return 0; } statement) through the program until the first line in which i found the seg fault occurring, and this is the line.
The whole thing is: The it works until the bold starts, then stops, i.e. the seg fault appears to occur at the last unbold line. There may be commented out areas in /* */ which should be included in the final code. Sorry for the italics, in a bit of a rush to change all the [ i ].
// Program to output the (xyz) origins and (lmn) directional cosines of the Cherenkov cone given off by a set particle
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <cassert>
using namespace std;
int main(){
//Initialise the myriad of variables & constants that will be needed
//Componemts for Sellmeier formula for Herasus SUPRASIL
//Constants
double A1 = 0.4731155913;
double A2 = 0.631039719;
double A3 = 0.987685322;
double B1 = 0.01995717;
double B2 = 0.0041280922;
double B3 = 98.7685322;
double PI = 3.14159265358979323846;
//Variables, lambda must be in microns
double minlambda = 0.0;
double maxlambda = 0.0;
int lambdasteps = 1;
double lambdas[30];
double refindex[30];
//Components for the mementa of the incoming particle
//Variables, momenta in GeV/c
double maxmomentum = 0.0;
double minmomentum = 0.0;
int momentumsteps = 1;
double momenta[30];
double betas[30];
//Components for the trajectory of the particle
double maxangle = 0.0;
double minangle = 0.0;
double maxanglein = 0.0;
double minanglein = 0.0;
int anglesteps = 1;
double angles[30];
int numberofpoints = 1;
double x[30];
double y[30];
double z[30];
//Cherenkov Angle
double cherenkovangles[30][30];
//Define the particles as far as I care, masses in GeV/c^2
double masschargedpion = 0.13957018;
double masschargedkaon = 0.493677;
double m0;
int mchoose = 0;
//Additional variables to be used
double radius[30][30];
double xprimevec[30][30][30][30];
double yprimevec[30][30][30][30];
double zprimevec[30][30][30][30];
double maxdirecangle=0.0;
double mindirecangle=0.0;
double maxdirecanglein=0.0;
double mindirecanglein=0.0;
int direcanglesteps=1;
double direcangles[30];
double vprime[3][30][30][30][30];
int xhat[3];
int yhat[3];
int zhat[3];
//Set up xhat,yhat and zhat
for(int i=0;i<=2;i++){
if(i==0){
xhat[i]=1;
}
else{
xhat[i]=0;
}
}
//Read in user defined variables
cout << "Please input the minimum wavelength, in microns, for use in the Sellmeier Equation: " << endl; cin >> minlambda;
cout << "And the similar maximum wavelength in microns: " << endl;
cin >> maxlambda;
cout << "And the number of steps to divide the range: " << endl;
cin >> lambdasteps;
cout << "Please maximum trajectory angle in degrees: " << endl;
cin >> maxanglein;
cout << "And the minimum angle in degrees: " <<endl;
cin >> minanglein;
cout << "And the number of steps to divide the range: " << endl;
cin >> anglesteps;
cout << "Please input the maximum momentum in GeV/c: " << endl;
cin >> maxmomentum;
cout << "And the minimum momentum in GeV/C: " << endl;
cin >> minmomentum;
cout << "And the number of steps to divide the range: " << endl;
cin >> momentumsteps;
cout << "Please input the minimum angular deviation of propagation from the +ve x-axis in degrees: " << endl;
cin >> mindirecangle;
cout << "And the maximum angle in degrees: " << endl;
cin >> maxdirecangle;
cout << "And the number of steps to divide the range: " << endl;
cin >> direcanglesteps;
cout << "Please press 0 for charged Pion and 1 for charged Kaon" << endl;
cin >> mchoose;
//STEP 1 ---- Use of Sellmeier formula to get Refractive index at given Lambdas
//Firstly calculate the wavelengths from given inputs
double lambdastepsize=((maxlambda-minlambda)/(lambdasteps));
//Secondly make a loop over the number of steps
for(int i=0; i<=(lambdasteps-1); i++){
//Third, fill an array of wavelengths
lambdas[i]=minlambda+((i)*lambdastepsize);
//Finally take those wavelengths and fire them into the Sellmeier equation to fill an array of refractive indices
refindex[i]=sqrt(1.0+((A1*(pow(lambdas[i],2.0)))/((pow(lambdas[i],2.0))-B1))+((A2*(pow(lambdas[i],2.0)))/((pow(lambdas[i],2.0))-B2))+((A3*(pow(lambdas[i],2.0)))/((pow(lambdas[i],2.0))-B3)));
}
//STEP 2 ---- Turn momenta into Cherenkov angle
//First calculate the momenta from given inputs
double momentastepsize=((maxmomentum-minmomentum)/(momentumsteps));
//loop over the number of steps
for(int j=0; j<=(momentumsteps-1);j++){
//fill array of momenta
momenta[j]=minmomentum+((j)*momentastepsize);
//calculate betas
betas[j]=sqrt((pow(momenta[j],2.0))/((pow(m0,2.0))+(pow(momenta[j],2.0))));
}
//Start a double nested loop to create the cherenkov angles involved
for(int i=0;i<=(momentumsteps-1);i++){
for(int k=0;k<=(lambdasteps-1);k++){
cherenkovangles[i][k]=acos((1.0)/(refindex[k]*betas[i]));
}
}
//STEP 3 ---- Calcualte trajectories of the particles
//Create array of angles of incidence for the incoming particles.
double anglestepsize=((maxangle-minangle)/(anglesteps));
double ystepsize=((10.000000-(-10.000000))/(numberofpoints));
for(int i=1;i<=anglesteps;i++){
angles[i]=minangle+((i-1)*anglestepsize);
}
//Create the arrays of the origins of the particles.
for(int i=1;i<=numberofpoints;i++){
x[i]=-25.000000;
y[i]=10.00000-((i-1)*ystepsize);
z[i]=tan(angles[i])*y[i];
}
//Set radius of unrotated circle
for(int i=1;i<=momentumsteps;i++){
for(int j=1;j<=lambdasteps;j++){
radius[i][j]=tan(cherenkovangles[i][j]);
}
}
//Fill the components of the matrices
double direcanglestepsize=((maxdirecangle-mindirecangle)/(direcanglesteps));
for(int a=1;a<=direcanglesteps;a++){
direcangles[a]=mindirecangle+((a-1)*direcanglestepsize);
for(int b=1;b<=anglesteps;b++){
for(int c=1;c<=momentumsteps;c++){
for(int d=1;d<=lambdasteps;d++){
xprimevec[a][b][c][d]=cos(direcangles[a])*tan(cherenkovangles[c][d]);
yprimevec[a][b][c][d]=-cos(angles[b])+sin(angles[b])*(sin(direcangles[a])*tan(cherenkovangles[c][d]));
zprimevec[a][b][c][d]=cos(angles[b])*(cos(direcangles[a])*tan(cherenkovangles[c][d])+sin(angles[b]));
}
}
}
}
//Putting the various data into a vector.
for(int a=1;a<=direcanglesteps;a++){
for(int b=1;b<=anglesteps;b++){
for(int c=1;c<=momentumsteps;c++){
for(int d=1;d<=lambdasteps;d++){
vprime[1][a][b][c][d]=xprimevec[a][b][c][d];
vprime[2][a][b][c][d]=yprimevec[a][b][c][d];
vprime[3][a][b][c][d]=zprimevec[a][b][c][d];
}
}
}
}
Have you considered the amount of memory used to allocate the 4 and 5 dimensional arrays of doubles? These are being allocated from the stack (vs. the heap). xprimevec[30][30][30][30] uses 30*30*30*30*sizeof(double) bytes of stack space. You have declared 6 of these 4 dimensional arrays and one 5 dimensional that is actually 3 additional 4-dimensional arrays of doubles. So, for a rough estimate, given sizeof(double) to be 8 and 9 arrays of 30^4 bytes, your code is trying to allocate well over 55Mb of stack space. I don't know if this is causing the segmentation fault but it is the first thing that jumps out at me.