3
\$\begingroup\$

I have written a code that performs numerical integrations—one using the Gauss method and the other using the trapezoidal method. For the Gauss method, I have five separate text files containing the weight-points (i.e., the zeros of the Laguerre polynomials) for different orders. These are necessary for the integration.

The struct defines my integrands, which depend on the temperature i and the chemical potential j. x is my integration variable and k the cases to calculate the integrands. Additionally, the function includes the necessary implementations for both the trapezoidal and Gauss integration methods.

In the main function, I use a loop to perform the integration for different temperatures and chemical potentials. The loop also contains expressions required for the calculations.

Can someone tell me why my program is taking so long to run, which operation is the slowest, and how I can efficiently improve it?

So far, I have optimized the program by introducing structs. However, it still takes a while to compute values.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <float.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// Constants

const double y = 2.0;
const double v = 1.0;
const double m = 938;
const double vorfac = 2.0 * M_PI * M_PI ;

// Gauss-Laguerre

const int n = 32; // Order of Gauss-Laguerre
double a = 0.0;     // Power of x, must be specified for each case in the function


typedef struct 
{
    double E , f_plus, f_minus, f_S_plus, f_S_minus ; 
} FuncBaseresult;



FuncBaseresult FuncBase(double x, double T, double mu)
{
    FuncBaseresult res;
    res.E = (sqrt(x*x + m*m ));
    double exp_E_T = exp(res.E/ T);
    double exp_mu_T = exp(mu / T);
    res.f_minus = 1.0 / (exp_E_T / exp_mu_T + v);
    res.f_plus = 1.0 / (exp_E_T * exp_mu_T + v);

    res.f_S_minus = - (mu - res.E) * exp( res.E / T + mu / T) * 1.0/T * 1.0/T * 1.0/(pow((v * exp(mu/T) + exp(res.E/T)) , 2)); 
     res.f_S_plus =  (mu + res.E) * exp( res.E / T + mu / T) * 1.0/T * 1.0/T * 1.0/(pow((  exp( res.E / T + mu / T) + v   ) , 2));

    return res;
}


double Function(int k, double x, double T, double mu) // Integrand
{
    FuncBaseresult base = FuncBase(x, T, mu);


    if (k == 0)
    {
        a = 2.0;
        return y/(vorfac) * x*x * (base.E-m) * (base.f_minus + base.f_plus); // Energy density 
    }

    if (k == 1)
    {
        a = 4.0;
        return  y /(vorfac) * 1.0 / 3.0 * x*x*x*x *  1.0/base.E *  (base.f_minus + base.f_plus); // Pressure
    }

    if (k == 2)
    {
        a = 4;
        return y /(vorfac) * 1.0 / 3.0 * x*x*x*x * 1.0/base.E * (base.f_S_minus + base.f_S_plus);
    }

    if (k == 3)
    {
        a = 2;
        return  y / (vorfac) * x*x * ( base.f_minus + base.f_plus); // Particle density
    }

    return 0;
}

double x_roots[32];  // Global speichern

void getroots() {
    char filename[10];
    sprintf(filename, "a=%d.txt", (int)a);
    
    FILE *file = fopen(filename, "r");
    if (file == NULL) {
        printf("Error opening file!\n");
        exit(1);
    }

    for (int i = 0; i < 32; i++) {
        fscanf(file, "%lf", &x_roots[i]);
    }
    fclose(file);
}


double fac(int n)
{
    double res = 1;
    for (int i = 1; i <= n; i++)
    {
        res *= i;
    }
    return res;
}

double laguerre(int n, double a, double x) {
    if (n == 0) return 1.0;
    if (n == 1) return 1.0 + a - x;

    double L0 = 1.0;
    double L1 = 1.0 + a - x;
    double L2 = 0.0;

    for (int i = 2; i <= n; i++) {
        L2 = ((2 * i - 1 + a - x) * L1 - (i - 1 + a) * L0) / i;
        L0 = L1;
        L1 = L2;
    }

    return L1;
}

double Gauss(int k, double *x, double T, double mu)
{
    double fac_n = fac(n);
    double fac_a = fac(a);
    double fac_na = fac(n+a);
    double res;
    res = 0;
    res *= Function(k, 0, T, mu);
    getroots(x);
    for (int i = 0; i < 32; i++)
    {
    double lag_val = laguerre(n + 1, a , x[i]);  // EINMAL berechnen
    res += Function(k, x[i], T, mu) * exp(x[i]) * x[i] * fac_na / 
          (fac_n * (n + 1) * (n + 1) * pow(lag_val, 2.0) * pow(x[i], a));
    }

    return res;
}

// Trapezoidal method

const double e = 1e-5; // Precision until stopping, method stops when new terms contribute less than this value

double d = 1e-3; // Interval size

double trapez(int k, double T, double mu)
{
    double res = 0;

    for (int i = 1; ;i++)
    {
        double f_val = Function(k, i * d, T, mu);

        if(fabs(f_val) <= e *res) break;
        res += d * f_val;
    }
    return res;
}

int main()
{
    FILE *file1 = fopen("test.txt", "w");
    if (file1 == NULL)
    {
        printf("Error opening the file!\n");
        exit(1);
    }

    double x[32];
    double temp0, temp1, temp2, temp3 ;
    double Euler;

    // j is mu and i is T

    for (double j = 1; j < 50000; j*=30) //mu
    {
        fprintf(file1, "%e ", j/m); // Start line with mu

        for (double i = 20; i <= 2000; i += 20) // Temperature
        {
            if (i == 20 || i == 1000  || i == 1500 )
            {
                if (j/i > 1e-5)
                {
                    temp0 = trapez(0, i, j);
                    temp1 = trapez(1, i, j);
                }
                else
                {
                    temp0 = Gauss(0, x, i, j);
                    temp1 = Gauss(1, x, i, j);
                }

                double k_f = sqrt(j*j - m*m); // j is mu
                double P_0;
                double e_0; 
                if (j < m)
                {
                    P_0 = 0; 
                    e_0 = 0;
                }
                else
                {
                    P_0 = y/(24.0*M_PI*M_PI) * (j * k_f * (k_f *k_f - 3.0/2.0 * m*m)  + 3.0/2.0 * pow(m ,4) * log( (k_f + j)/m ));
                    e_0 =  y/(2.0*M_PI*M_PI) * 1.0/24.0 * (-3*asinh(k_f/m) * pow(m,4)  + j * (3 * k_f * m*m + 6*pow(k_f , 3))  - 8 * pow(k_f , 3) *m );
                }

                if(temp0  <  1e-148 || temp1  < 1e-148 )
                {
                    temp0 =  8.170380e-170;
                    temp1 =  1.630096e-169 ; 
                }
                
                double gamma2 = ((temp1 - P_0)/(temp0 - e_0  )) + 1 ; 

                fprintf(file1, " T/m = %e gamma2 = %e ||",i/m , gamma2);
            }
        }

        fprintf(file1, "\n");
    }

    fclose(file1);
}

a=0.txt


0.04448936583326701388
0.23452610951961938568
0.57688462930187434097
1.07244875381808424386
1.72240877643981793454
2.52833670644600161381
3.49221327314977170175
4.61645676498042956126
5.90395853784575752599
7.35812653405240890692
8.98294165535760669172
10.78301623367992156943
12.76371176580493127517
14.93115838768134118197
17.29244179823459504064
19.85579651978946458257
22.63255043922785247901
25.63028835810287731078
28.86407464765523656069
32.35453780640911247701
36.10748356744834808296
40.18143213593468487943
44.42889069562876613873
49.27734571508544547669
54.27154725396455603459
59.94403361568824095684
65.93755634763675743670
72.70366493025684917484
80.18138178742300681279
88.73656778030174052674
98.82940647607993867041
111.75140449427495070722

a=1.txt



0.11125830781785717161
0.37321887563579131397
0.78564593354865541741
1.34945557597258414262
2.06595254242072190465
2.93682515070008509639
3.96416490359780615549
5.15049471449633600173
6.49880480363346535455
8.01259753882315450824
9.69594322587912316180
11.55354602478912795505
13.59083274319271872344
15.81411321911575207366
18.23009213725826782593
20.84849310402791644492
23.67515123618571948327
26.73438742869342021891
29.99552028898686018010
33.57676206461579226925
37.34146675218031674603
41.46495585451080501116
45.88785208864274522966
50.61186477968586672205
55.84742703528508656063
61.37005353675105112643
67.57416441994908495872
74.30608614857651161856
81.87696228978541057586
90.48403294680814212825
100.64522034615339407537
113.64498858551716864440

a=2.txt

0.19694392214666700536
0.52948786605016606721
1.01026981913835722793
1.64061619167524574969
2.42200673329678828338
3.35625823788337251941
4.44557318796410427808
5.69257571438529819119
7.10035054452460823171
8.67248886206972890989
10.41314853479853397289
12.32708584885750902060
14.41998155470366604902
16.69734167676762481847
19.16999706336005004914
21.83522688699483182972
24.73213241077598212314
27.78526342547449345943
31.19219783012959368307
34.70691588805129157436
38.66256189358491468511
42.80051011906167701682
47.08750526143398928980
52.18367944128409874338
57.16929227843460381564
63.00419269308171266175
69.05512238051676376926
75.96022652445722656012
83.54518851546569635502
92.22284353379730248434
102.44781654325494457680
115.52491081961949248580

a=3.txt

0.29961872904924113925
0.70198106535399051875
1.24974569814540892310
1.94514382444181599396
2.78994869155615976553
3.78614305436027143159
4.93605294225746238368
6.24240683549429142829
7.70838397617374493365
9.33766076211349371761
11.13449041432764907711
13.10369376236954686021
15.25092379849237644862
17.58290560002000546547
20.10295780865619619249
22.84125955809728836243
25.74399385028829456701
32.24991591099860244185
36.02127473917757782829
43.96048256063561865403
48.61417396672651847211
58.78464627570348710606
64.41104405231406815346
70.68058601977669752614
77.53803838832618566812
85.21629249041217235572
93.94733619882016739666
104.23846159852307380334
117.39174740382649986259
39.76478080625362565570
53.26931493053049848641
28.98039246758906983814

a=4.txt

0.41790923492851933529
0.88953757977417469149
1.50313435358155333965
2.26227724680885167174
3.16916090796869021062
4.22598142404881293999
5.43520757619423022788
6.79968010050789395393
8.32267392651077386745
10.00795498422072071776
11.85983746177703679336
13.88318378640620309739
16.08389351549142176623
18.46814636381402863208
23.81677094107534387035
26.81166653707526492667
30.06799305570848446223
33.37113780215095459880
40.76480747240096746964
45.63178176901075744354
72.06876424060635599744
119.24607803365869074241
21.04413848219094163028
66.17739535280246343518
37.32268320943082784424
55.39666901945717114586
95.66257883144125173658
49.44734140119586385254
59.78453975634411676765
86.86244827669744950072
79.18793460362063285629
106.01725817767682258363

I have now tried to implement the suggestions and optimized my code. The code is running faster, but some integrands are being calculated incorrectly, for example when k==3. Did I make any mistakes in the implementation?

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <float.h>
#include <omp.h>  // OpenMP für Parallelisierung

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

// Konstanten

const double y = 2.0;
const double v = 1.0;
const double m = 930;
const double vorfac = 2.0 * M_PI * M_PI;

const int n = 32; // Gauss-Laguerre Ordnung
double a = 0.0;

double x_roots[32];  // Save weights
double fac_lookup[33]; // save fac!

typedef struct {
    double E, f_plus, f_minus, f_S_plus, f_S_minus;
} FuncBaseresult;


// save fac
void init_fac_table() {
    fac_lookup[0] = 1.0;
    for (int i = 1; i <= 32; i++) {
        fac_lookup[i] = fac_lookup[i - 1] * i;
    }
}


void getroots(double a_value) { 
    char filename[10];
    sprintf(filename, "a=%d.txt", a_value);  
    
    FILE *file = fopen(filename, "r");
    if (file == NULL) {
        printf("Error opening file: %s\n", filename);
        exit(1);
    }

    for (int i = 0; i < 32; i++) {
        fscanf(file, "%lf", &x_roots[i]);
    }
    fclose(file);
}

FuncBaseresult FuncBase(double x, double T, double mu) {
    FuncBaseresult res;
    res.E = sqrt(x*x + m*m);
    
    double exp_E_T = exp(res.E / T);
    double exp_mu_T = exp(mu / T);
    double exp_sum = exp_E_T * exp_mu_T;
    double exp_diff = exp_E_T / exp_mu_T;

    res.f_minus = 1.0 / (exp_diff + v);
    res.f_plus = 1.0 / (exp_sum + v);

    double denom1 = (v * exp_mu_T + exp_E_T);
    double denom2 = (exp_sum + v);


    
    res.f_S_minus = - (mu - res.E) * exp(exp_E_T + exp_mu_T) * (1.0 / (T*T)) / (denom1 * denom1);
    res.f_S_plus  =  (mu + res.E) * exp(exp_E_T + exp_mu_T) * (1.0 / (T*T)) / (denom2 * denom2);

    return res;
}

double Function(int k, double x, double T, double mu) {
    FuncBaseresult base = FuncBase(x, T, mu);

    if (k == 0) {
        a = 2.0;
        return y / vorfac * x*x * (base.E - m) * (base.f_minus + base.f_plus);
    }
    if (k == 1) {
        a = 4.0;
        return y / (vorfac * 3.0) * x*x*x*x / base.E * (base.f_minus + base.f_plus);
    }
    if (k == 2) {
        a = 4.0;
        return y / (vorfac * 3.0) * x*x*x*x / base.E * (base.f_S_minus + base.f_S_plus);
    }
    if (k == 3) {
        a = 2.0;
        return y / vorfac * x*x * (base.f_minus + base.f_plus);
    }

}



// Laguerre-recursive
double laguerre(int n, double a, double x) {
    double L0 = 1.0, L1 = 1.0 + a - x, L2;
    for (int i = 2; i <= n; i++) {
        L2 = ((2 * i - 1 + a - x) * L1 - (i - 1 + a) * L0) / i;
        L0 = L1;
        L1 = L2;
    }
    return L1;
}

// Gauss-Quad
double Gauss(int k, double T, double mu) {
    double fac_n = fac_lookup[n];
    double fac_a = fac_lookup[(int)a];
    double fac_na = fac_lookup[n + (int)a];

    double res = 0;
    for (int i = 0; i < 32; i++) {
        double lag_val = laguerre(n + 1, a, x_roots[i]);
        res += Function(k, x_roots[i], T, mu) * exp(x_roots[i]) * x_roots[i] * fac_na / 
               (fac_n * (n + 1) * (n + 1) * pow(lag_val, 2.0) * pow(x_roots[i], a));
    }
    return res;
}


// Trapez
const double e = 1e-5;
double d = 1e-3;

double trapez(int k, double T, double mu) {
    double res = 0, prev_res = 0;
    double h = d;

    for (int i = 1; ; i++) {
        double f_val = Function(k, i * h, T, mu);
        res += h * f_val;

        if (fabs(res - prev_res) < e * res) break;
        prev_res = res;
        h *= 0.9;
    }
    return res;
}





int main() {
    FILE *file1 = fopen("test.txt", "w");
    if (!file1) {
        printf("Error opening the file!\n");
        return 1;
    }

    init_fac_table();
    getroots(0); 

    double temp0;
    double temp1;
    #pragma omp parallel for schedule(dynamic) // Parallelisierung
    for (double j = 1; j < 1e6; j *= 10) {  // mu
    fprintf(file1, "%e ", j ); // Start line with mu

    for (double i = 20; i <= 2000; i += 20) { // Temperature
        if (i == 20 || i == 400  || i == 1000  || i == 2000){

        double k_f = sqrt(j * j - m * m);
        double P_0 = (j >= m) ? y / (24.0 * M_PI * M_PI) * (j * k_f * (k_f * k_f - 3.0 / 2.0 * m * m) + 3.0 / 2.0 * pow(m, 4) * log((k_f + j) / m)) : 0 ; 
        double e_0 = (j>=m) ? y / (2.0 * M_PI * M_PI) * 1.0 / 24.0 * (-3 * asinh(k_f / m) * pow(m, 4) + j * (3 * k_f * m * m + 6 * pow(k_f, 3)) - 8 * pow(k_f, 3) * m) : 0 ;
        
        
        double P_norm = 2 * (7.0 / 4.0 * M_PI * M_PI / 90 * pow(i, 4) + 1.0 / 12.0 * pow(j * i, 2) + 1.0 / 24.0 * 1 / pow(M_PI, 2) * pow(j, 4)); //Faktor 2 am Anfang, wegen Normierung von Teilchen-Antiteilchen

         double S_norm = 2*(7.0 * M_PI * M_PI / 90 * pow(i, 3) + 1.0 / 6.0 * pow(j, 2) * i);

        double n_norm;
        if (j == 0)
        {
            n_norm = 2.0;
        }
        else
        {
            n_norm = 2*(1.0 / 6.0 * j * pow(i, 2) + 1.0 / 6.0 * 1 / (pow(M_PI, 2))*pow(j, 3));
        }
        

        double temp0 = (j / i > 1e-5) ? trapez(0, i, j) : Gauss(0, i, j);  
        double temp1 = (j / i > 1e-5) ? trapez(1, i, j) : Gauss(1, i, j); 
        double temp3 =  (j / i > 1e-5) ? trapez(3, i, j) : Gauss(3, i, j);
        temp3=temp3/n_norm; 
        double gamma2 = (temp1 - P_0) / (temp0 - e_0) + 1;  

        fprintf(file1, " n = %e  ||", temp3);
    }
    }

    fprintf(file1, "\n");
}


    fclose(file1);


    return 0;

}
```
\$\endgroup\$
6
  • 1
    \$\begingroup\$ What did profiling the program tell you about where the time went? \$\endgroup\$ Commented Mar 22 at 17:06
  • 3
    \$\begingroup\$ Something is wrong with your code. 1 getroots reads the a*.txt files into x_roots, which is never referenced. It is as good as not read them at all. 2 I tried to run your code, and found that trapez never reaches break. Back to the drawing board I guess. VTC. \$\endgroup\$ Commented Mar 22 at 22:58
  • \$\begingroup\$ The code really never worked correctly and that made the question off-topic for the Code Review site. \$\endgroup\$ Commented Mar 23 at 21:40
  • \$\begingroup\$ @pacmaninbw The close reason says that code must work correctly to to the best of the authors knowledge. So then should questions that do not work correctly due to a subtle bug that the author did not know about still be closed? \$\endgroup\$ Commented Mar 23 at 22:52
  • \$\begingroup\$ @CPlus The original code used an uninitialized array of double as the roots, it isn't clear how it ever caculated anything. I did not close the question, it was closed by the community. I do my best not to answer off-topic questions. \$\endgroup\$ Commented Mar 23 at 23:18

2 Answers 2

Reset to default
3
\$\begingroup\$

Some ideas to improve performance:

If you can get by with float, use float objects, constants and float functions like expf().

Use a code profiler to identify hot spots.

First look to the function calls and save results. A compiler may not see this optimization and so benefits with a little help.

res.f_S_minus = - (mu - res.E) * exp( res.E / T + mu / T) * 1.0/T * 1.0/T * 
    1.0/(pow((v * exp(mu/T) + exp(res.E/T)) , 2)); 
res.f_S_plus  =   (mu + res.E) * exp( res.E / T + mu / T) * 1.0/T * 1.0/T * 
    1.0/(pow((  exp( res.E / T + mu / T) + v   ) , 2));

may improve with:

double ex = exp( res.E / T + mu / T);
res.f_S_minus = - (mu - res.E) * ex * 1.0/T * 1.0/T * 1.0/(pow((v * exp(mu/T) + exp(res.E/T)) , 2)); 
res.f_S_plus  =   (mu + res.E) * ex * 1.0/T * 1.0/T * 1.0/(pow((  ex + v   ) , 2));

Consider saving 1.0/T * 1.0/T also for similar reasons.

// pow(k_f , 3)

may improve with:

(k_f * k_f * k_f)

// pow(lag_val, 2.0)
(lag_val*lag_val)

// Save `m*m` from before.
double mm = m*m;

// and then
// pow(m,4)
mm*mm

2nd, look to trading division for multiplication. (I did not find any.)

Review candidate operations to reduce errors.

When j is near m, performing the subtraction before multiplying is more accurate. It also trades 1 product, 2 sums for 2 products, 1 sum.

// sqrt(j*j - m*m);
sqrt((j-m) * (j+m));

Note that trading sqrt(x*x + m*m ); for hypot(x, m); is more accurate, but can be slower.

Replace iterative code with table lookup as n is small for fac() (32 + 4).

#include <assert.h>

double fac(unsigned n) {
  static const double f[] = {
      1, //
      1, //
      1.0*2, //
      1.0*2*3, //
      1.0*2*3*4, //
      1.0*2*3*4*5, //
      1.0*2*3*4*5*6, // Let the compile do the math at compile time.
      // ...
      };
  assert(n < sizeof f/sizeof f[0]);
  return f[n];
}

If you do not want to load a table up ahead of time, consider saving the results to avoid re-computing.

Replace double a with unsigned a.

\$\endgroup\$
2
  • \$\begingroup\$ float computations are not necessarily faster than double computations. Intel chips for example always compute internally using long double for precision, unless you’re using SIMD instructions. \$\endgroup\$ Commented Mar 23 at 13:17
  • \$\begingroup\$ I agree that your fac() function is better, but the getroots() function which opens a file and reads from a file is called from Gauss() which us called from a nested loop in main(). The single most effective optimization would be to call getroots() once in main(). \$\endgroup\$ Commented Mar 23 at 15:29
4
\$\begingroup\$

General Observations

If the function getroots() can't open the input file the program does not execute any clean up code.

Optimization

An Optimization that would definitely speed up the program is to call the getroots() function only once and pass the input array into the Gauss() function as necessary. File I/O is very expensive in terms of execution time.

When helping to optimize code it helps us to know what operating system, processor and compiler you are using. For instance are you building and running this on Linux or on Windows? Are you using the GNU compilers, Clang or Microsoft Visual Studio? What version of C are you using? If you are using Linux and GCC are you compiling -O3 which is the most optimized generated code. Debuggable code is generally slower.

It also helps if you can provide profiling information. On Linux using gcc compile the program using gcc -g -pg -o myprogram myprogram.c.

  • -g: Includes debugging information, which is useful for associating profiling data with source code lines.
  • -pg: Enables profiling instrumentation, which adds code to collect execution statistics.

Run the program. This generates a gmon.out file containing the profiling data. Analyze the profiling data using gprof.

Run gprof gprof myprogram gmon.out.

This command produces a textual report with information about:

Call graph:
Shows the function call hierarchy and the number of times each function was called.

Flat profile:
Lists functions and the time spent in each, as well as the percentage of total execution time.

Let the compiler help you write better code.

If you are compiling using the GCC compiler use the compiler options -Wall -Wextra and -Werror. When compiling this I found that there are multiple variables that are not used in the program:

In the function Gauss()

    double fac_a = fac(a);

In main() there are 3 unused variables, temp2, temp3 and Euler.

Avoid Global Variables

It is very difficult to read, write, debug and maintain programs that use global variables. Global variables can be modified by any function within the program and therefore require each function to be examined before making changes in the code. In C and C++ global variables impact the namespace and they can cause linking errors if they are defined in multiple files. The answers in this stackoverflow question provide a fuller explanation.

Magic Numbers

While many constants are already defined, there are more magic numbers that really should be turned into constants. In main():

    for (double j = 1; j < 50000; j*=30) //mu
    {
        fprintf(file1, "%e ", j/m); // Start line with mu

        for (double i = 20; i <= 2000; i += 20) // Temperature
        {
            if (i == 20 || i == 1000  || i == 1500 )

The numbers 50000, 30, 20, 2000, 1000 and 1500 are all magic number and a named constant would make it easier to undstand and maintain the code.

Variables and Function Names

Rather than have a comment that explains j is mu and i is temperature, the variables should be named mu and temperature.

The function name Gauss() is descriptive, the function getroots() is descriptive, the function name trapez() is somewhat descriptive I might have used the full word, the same for the function fac(). The functions FuncBase() and Function() are not descriptive at all.

It isn't clear what the variable k in either Gauss() or trapez().

Use Code That Can be Easily Expanded

In the function Function() I would replace all the if statements with a switch/case statement.

double Function(int k, double x, double T, double mu) // Integrand
{
    FuncBaseresult base = FuncBase(x, T, mu);

    switch(k)
    {
        default:
        {
            fprintf(stderr, "Unknown Option in Function\n.");
            return 0;
        }
        case 0:
        {
            a = 2.0;
            return y/(vorfac) * x*x * (base.E-m) * (base.f_minus + base.f_plus); // Energy density 
        }
        case 1:
        {
            a = 4.0;
            return  y /(vorfac) * 1.0 / 3.0 * x*x*x*x *  1.0/base.E *  (base.f_minus + base.f_plus); // Pressure
        }
        case 2:
        {
            a = 4;
            return y /(vorfac) * 1.0 / 3.0 * x*x*x*x * 1.0/base.E * (base.f_S_minus + base.f_S_plus);
        }
        case 3:
        {
            a = 2;
            return  y / (vorfac) * x*x * ( base.f_minus + base.f_plus); // Particle density
        }
    };
}

This is another place where there are magic numbers.

Complexity

The function main() is too complex (does too much). As programs grow in size the use of main() should be limited to calling functions that parse the command line, calling functions that set up for processing, calling functions that execute the desired function of the program, and calling functions to clean up after the main portion of the program.

There is also a programming principle called the Single Responsibility Principle that applies here. The Single Responsibility Principle states:

that every module, class, or function should have responsibility over a single part of the functionality provided by the software, and that responsibility should be entirely encapsulated by that module, class or function.

\$\endgroup\$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.