Fixing wrong radius values on table output?

Blueshark1 (29)
Hello,
For a project I am supposed to have my code output a table as temperature and radius change going from the core of a star to its surface. So far the code can run, and it shows the change in temperature from the core to the surface with the given inputs, but it does not show a proper change in radius, as the temperature goes down, the radius should go up, but my program only outputs -inf and zeroes for the radii, or as the temperature drops, so does the radius, which doesn't make sense.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
#include<iostream>		//Required for cin, cout
#include<cmath>			
using namespace std;
int
main ()
{
  double G(6.67*pow(10,-11)), rad,k(1.38*pow(10, -23)), RadCore, meanmolweight(0.6060), centdens(160), protonmass(1.67*pow(10,-27)), x_incr, temp, estmass, a, b, new_x;	// Declare objects in main()
  int Number, loop;
  cout <<
    "Enter the radius of the star in meters\n";
  cin >> rad;  
  cout <<
    "Enter the radius of the core in meters\n";
    cin >> RadCore;
  cout <<
    "Enter the estimated mass of the star in kg\n";
  cin >> estmass; 

  cout << "Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted \n";
  cin >> a >> b;
  cout <<
    "Enter the number of points you want over the range of 'a' to 'b'\n";
  cin >> Number;
  x_incr = (b - a) / Number;	// increment between the points based on number wanted
  loop = int (b - a) / x_incr;	// loop gives the for loop repetitions equal to number points wanted
  cout << loop;			// check that shows loop and Number are the same
  cout.setf (ios::fixed);
  cout.precision (2);		// Set Formats
  cout << "Temperature and Radius (new_x) \n";
  for (int k = 0; k <= loop; k++)
    {
      new_x = a + k * x_incr;
      temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*pow(10,2)*(1+(RadCore/rad)-((19*RadCore*pow(10,2))/(5*rad*pow(10,2)))+(((9*RadCore*pow(10,3))/(5*rad*pow(10,3)))));
      cout << new_x << " " << temp << endl;	//echo vales to screen
    }
  return 0;	}		// Exit program.
Last edited on
salem c (3685)
Fix the code tags please!
Peter87 (11178)
I see that you divide by (36*k) but if k is equal to zero that means you'll be dividing by zero which would result in an "infinite" floating-point value.
Last edited on
lastchance (6980)
@Blueshark1,
At the start of your code you use k for what is presumably the Boltzmann constant:
k(1.38*pow(10, -23))


Then, unfortunately, you also use it for a loop variable:
1
2
3
4
5
 
for (int k = 0; k <= loop; k++)
{
   new_x = a + k * x_incr;
// ...  
}


The latter will hide the (intended) value of the Boltzmann constant and lead, as @Peter87 says, to potential division by 0.


Please:
- PUT YOUR CODE IN CODE TAGS!!!
- dont't write things like k(1.38*pow(10, -23)) Just use scientific notation k{1.38e-23}
- don't state the bleedin' obvious! e.g. "// loop gives the for loop repetitions equal to number points wanted" or "// Declare objects in main()"
- use a decent value for PI, not 3.14
- give us some sample input.
Last edited on
seeplus (6457)
To use code tags:


[code]
// Paste formatted code here
[/code]


which gives:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
 
#include <iostream> //Required for cin, cout
#include <cmath>

using namespace std;

int main() {
	double G(6.67 * pow(10, -11)), rad, k(1.38 * pow(10, -23)), RadCore, meanmolweight(0.6060), centdens(160), protonmass(1.67 * pow(10, -27)), x_incr, temp, estmass, a, b, new_x; // Declare objects in main()
	int Number, loop;

	cout << "Enter the radius of the star in meters\n";
	cin >> rad;
	cout << "Enter the radius of the core in meters\n";
	cin >> RadCore;
	cout << "Enter the estimated mass of the star in kg\n";
	cin >> estmass;
	cout << "Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted \n";
	cin >> a >> b;
	cout << "Enter the number of points you want over the range of 'a' to 'b'\n";
	cin >> Number;

	x_incr = (b - a) / Number; // increment between the points based on number wanted
	loop = int(b - a) / x_incr; // loop gives the for loop repetitions equal to number points wanted

	cout << loop; // check that shows loop and Number are the same
	cout.setf(ios::fixed);
	cout.precision(2); // Set Formats
	cout << "Temperature and Radius (new_x) \n";

	for (int k = 0; k <= loop; k++) {
		new_x = a + k * x_incr;
		temp = ((5 * 3.14 * G * estmass * protonmass) / (36 * k)) * centdens * rad * pow(10, 2) * (1 + (RadCore / rad) - ((19 * RadCore * pow(10, 2)) / (5 * rad * pow(10, 2))) + (((9 * RadCore * pow(10, 3)) / (5 * rad * pow(10, 3)))));
		cout << new_x << " " << temp << endl; //echo vales to screen
	}
} // Exit program.


Also note that pow(10, 3) is 1000 and pow(10, 2) is 100.

cmath has pi already defined as a constant M_PI (you might need #define _USE_MATH_DEFINES before #include <cmath>).

Also note that if you are using at least C++20, then pi is defined in header <numbers>
https://en.cppreference.com/w/cpp/numeric/constants
Last edited on
Blueshark1 (29)
Here is my updated code, I was using the pow function incorrectly, but now all of my radii are the same. I apologize for my lack of knowledge, I haven't used C++ in a long time.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
 
#include<iostream>		//Required for cin, cout
#include<cmath>			
using namespace std;
int
main ()
{
  double G(pow(6.67,-11)), rad,kb(pow(1.38,-23)), RadCore, meanmolweight(0.6060), centdens(160), protonmass(pow(1.67,-27)), x_incr, temp, estmass, a, b, new_x;	// Declare objects in main()
  int Number, loop;
  cout <<
    "Enter the radius of the star in km\n";
  cin >> rad;  
  cout <<
    "Enter the radius of the core in km\n";
    cin >> RadCore;
  cout <<
    "Enter the estimated mass of the star in kg\n";
  cin >> estmass; 

  cout << "Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted \n";
  cin >> a >> b;
  cout <<
    "Enter the number of points you want over the range of 'a' to 'b'\n";
  cin >> Number;
  x_incr = (b - a) / Number;	// increment between the points based on number wanted
  loop = int (b - a) / x_incr;	// loop gives the for loop repetitions equal to number points wanted
  cout << loop;			// check that shows loop and Number are the same
  cout.setf (ios::fixed);
  cout.precision (2);		// Set Formats
  cout << "Temperature and Radius (new_x) \n";
  for (int k = 0; k <= loop; k++)
    {
      new_x = a + k * x_incr;
      temp = ((5*M_PI*G*estmass*protonmass)/(36*kb))*centdens*pow(rad,2)*(1+(RadCore/rad)-((19*pow(RadCore,2))/(5*pow(rad,2)))+(((9*pow(RadCore,3))/(5*pow(rad,3)))));
      cout << new_x << " " << temp << endl;	//echo vales to screen
    }
  return 0;	}		// Exit program.



/*Enter the radius of the star in km
70000000
Enter the radius of the core in km
10000 
Enter the estimated mass of the star in kg
1.9e+20
Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted 
20000 
6000
Enter the number of points you want over the range of 'a' to 'b'
10 
10Temperature and Radius (new_x) 
20000.00 89420920672767955991265280.00
18600.00 89420920672767955991265280.00
17200.00 89420920672767955991265280.00
15800.00 89420920672767955991265280.00
14400.00 89420920672767955991265280.00
13000.00 89420920672767955991265280.00
11600.00 89420920672767955991265280.00
10200.00 89420920672767955991265280.00
8800.00 89420920672767955991265280.00
7400.00 89420920672767955991265280.00
6000.00 89420920672767955991265280.00 */
lastchance (6980)
Well of course they come out the same: there is nothing in temp that changes with the iteration loop variable or new_x. Did you mean to replace rad by new_x or something?

BTW, I'd be very careful how you mix up units: the implied length units in both G and k are metres, but all your radii are being entered in km.
Last edited on
dhayden (5795)
Following on what seeplus said, 6.67 * pow(10,-11) can be written as 6.67E-11

Your heading is
cout << "Temperature and Radius (new_x) \n";

and print
cout << new_x << " " << temp << endl; //echo vales to screen

Do you see the problem here?

Your temperature calculation doesn't depend on the new_x, so it's the same for every iteration through the loop. It can also be greatly simplified:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*pow(10,2)*(1+(RadCore/rad)-((19*RadCore*pow(10,2))/(5*rad*pow(10,2)))+(((9*RadCore*pow(10,3))/(5*rad*pow(10,3)))));

temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*100*(1+(RadCore/rad)-((19*RadCore*100)/(5*rad*100))+(((9*RadCore*1000)/(5*rad*1000))));

temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*100*(1+(RadCore/rad)-
((1900*RadCore)/(500*rad))+(((9000*RadCore)/(5000*rad))));

temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*100*
       (1+(RadCore/rad) - ((19./5)*(RadCore/rad) + (9./5)*(RadCore/rad));

temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*100*
       (1 + (RadCore/rad) * (1 - 19./5 + 9./5));

temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*100*
       (1 + (RadCore/rad) * (-1));

temp = ((5*3.14*G*estmass*protonmass)/(36*k))*centdens*rad*100*
       (1 - (RadCore/rad));

Is that right, or is your formula wrong? Poking around, I think I found the formula you're using in equation 34 at [url]https://www.astro.utu.fi/~cflynn/Stars/l3.html[/url]. In that case, your code is wrong. Perhaps do this:

1
2
3
4
5
6
7
8
9
 
constexpr double PI = 3.141592653589793;  // Don't use 3.14 unless you want the accuracy of a slide rule.
...
   for (int k = 0; k <= loop; k++) {
        new_x = a + k * x_incr;
        double fr = new_x / rad;  r/R in the equation
        temp = (5*PI*G*protonmass*meanmolweight / (36*K) * centdens * rad * rad *
                (1 + fr - (19./5.)*fr*fr + (9./5.)*fr*fr*fr));
        cout << new_x << '\t' << temp << '\n';   //echo vales to screen
    }

Blueshark1 (29)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
 
#include<iostream>		//Required for cin, cout
#include<cmath>			
using namespace std;
int
main ()
{
  double G(pow(6.67,-8)), rad,kb(pow(1.38,-16)), RadCore, meanmolweight(0.6), protonmass(pow(1.67,-24)), x_incr, temp, estmass, a, b, new_x;	// Declare objects in main()
  int Number, loop;
  cout <<
    "Enter the radius of the star in cm\n";
  cin >> rad;  
  cout <<
    "Enter the radius of the core in cm\n";
    cin >> RadCore;
  cout << "Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted \n";
  cin >> a >> b;
  cout <<
    "Enter the number of points you want over the range of 'a' to 'b'\n";
  cin >> Number;
  x_incr = (b - a) / Number;	// increment between the points based on number wanted
  loop = int (b - a) / x_incr;	// loop gives the for loop repetitions equal to number points wanted
  cout << loop;			// check that shows loop and Number are the same
  cout.setf (ios::fixed);
  cout.precision (2);		// Set Formats
  cout << "Temperature and Radius (new_x) \n";
    for (int k = 0; k <= loop; k++) 

    { 

      rad = a + k * x_incr; 

      temp = ((5*M_PI*G*meanmolweight*protonmass)/(36*kb))*((RadCore-rad)/(rad*RadCore))*pow(rad,2)*(1+(RadCore/rad)-((19*pow(RadCore,2))/(5*pow(rad,2)))+(((9*pow(RadCore,3))/(5*pow(rad,3))))); 

      cout << rad << " " << temp << endl;	//echo vales to screen 

    } 
  return 0;	}		// Exit program.




/*Enter the radius of the star in cm
7e+10
Enter the radius of the core in cm
7e+9
Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted 
1.5e+7
6000
Enter the number of points you want over the range of 'a' to 'b'
10
10Temperature and Radius (new_x) 
15000000.00 142336.57
13500600.00 175826.01
12001200.00 222653.77
10501800.00 290965.83
9002400.00 396225.88
7503000.00 570794.07
6003600.00 892104.33
4504200.00 1585964.76
3004800.00 3566049.39
1505400.00 14216877.34
6000.00 895560636743.35 */


Thank you for your response, this is my newer code after changing some things around with the units and such, the code works, but only for a certain range then the formula no longer works (having the code output around 100 points, the orders of magnitude shoot up), if I only wanted the data table to output the range from 60,000K to 200,000K for example, what would I change? But I will try to change what I currently have with your advice.
Last edited on
jonnin (11331)
you just change the loop to match what you want.
probably just put a condition on the cout area:
if(something >= 60000 && something <= 200000)
{
cout as you have it
if( something > 200000) break; //maybe this to stop looping?
}

you can do it even cleaner if you back-solve for the start and end loop values that yield the desired values, assuming your function is increasing (I think it would be) and easy to work with? Then it would just be for (.. start to end with some increment) but the above is a simple approach that is probably good enough.
TheIdeasMan (6781)
....I was using the pow function incorrectly,...


As stated twice now by others, we don't want you to use the pow function at all. pow(6.67,-8) is 6.67-8 not 6.67 * 10-8 Just use G {6.67e-8} You seem to manage to do this in your inputs.

The use of the pow function for squaring is inefficient, just multiply them. When you do, store them as a separate variables, it will save you repeating yourself in the long expression, and save the compiler evaluating the same thing over and over.

((5*M_PI*G*meanmolweight*protonmass)/(36*kb))
That is a constant value, calculate it first, store it in a variable then put it in the larger expression. With the 36, I always write that as 36.0 to emphasise that it is a double as all the other values; it saves the compiler from doing an implicit cast to double from int. Doing that can help stop annoying bugs, such as integer division.

"Enter the radius of the star in cm\n";

I am intrigued, what type of nuttiness is this: who use centimetres in astrophysics? I see you have changed it from metres to km to cm, but ignored lastchance's advice about using metres. Are you messing with us?

Last edited on
seeplus (6457)
Reformatted (NOTE I don't know if the equations or the constants are correct!):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
 
#include <iostream>
#include <cmath>
using namespace std;

int main() {
	const double G { 6.6743e-8 }, kb { 1.38e-16 }, meanmolweight { 0.6 }, protonmass { 1.67262e-24 };
	double a {}, b {}, rad {}, RadCore {};
	int Number {};

	cout << "Enter the radius of the star in cm\n";
	cin >> rad;

	cout <<	"Enter the radius of the core in cm\n";
	cin >> RadCore;

	cout << "Enter endpoint range a and b (a<b) (First the core and then surface temperature of the star) of the function to be plotted \n";
	cin >> a >> b;

	cout << "Enter the number of points you want over the range of 'a' to 'b'\n";
	cin >> Number;

	const auto x_incr { (b - a) / Number };			// increment between the points based on number wanted
	const auto loop { int((b - a) / x_incr) };		// loop gives the for loop repetitions equal to number points wanted

	cout << loop << '\n';					// check that shows loop and Number are the same

	cout.setf(ios::fixed);
	cout.precision(2);		// Set Formats
	cout << "Temperature and Radius (new_x) \n";

	const auto calc { (5 * M_PI * G * meanmolweight * protonmass) / (36 * kb) };
	const auto radCore2 { RadCore * RadCore };
	const auto radCore3 { radCore2 * RadCore };

	for (int k = 0; k <= loop; k++) {
		rad = a + k * x_incr;

		const auto rad2 { rad * rad };
		const auto rad3 { rad2 * rad };

		const auto temp { calc * ((RadCore - rad) / (rad * RadCore)) * rad2 * (1 + (RadCore / rad) - ((19 * radCore2) / (5 * rad2)) + (((9 * radCore3) / (5 * rad3)))) };

		cout << temp << " " << rad << '\n';
	}
}


It can be seen that rad is defined L7 and value obtained L11. However this value isn't used. rad is re-calculated L36 and this value is used on subsequent lines.

Should be rad variable on L36 be called a different name? Do you need 2 variables - one input and one calculated? Which is used where in the calculation L41?

Also, shouldn't rad and temp be reversed on L43 to match the heading on L29?

L16 - you ask for a and b input with a < b. But with the sample input a > b ??
Last edited on
Blueshark1 (29)
The rad variable used in the temperature equation in the for loop should use the rad value from the input. I did try changing the variables, but it ended up messing up the output even more. As for why it is in cm, I genuinely don't know. The equation I am using from the link mentioned above has the units in cm and grams so I kept it as is. I did simplify the equation a bit by calculating the constants in the first part (like shown in seeplus"s code). However, I want to only have it output the results between the temperatures 200,000K and 60,000K. How would I apply an if statement into my for loop for that output only?
seeplus (6457)
Assuming that the displayed temp is in K, then in my code above insert before L43:

 
 
if (temp >= 60'000 && temp <= 200'000)


which gives you:

1
2
 
if (temp >= 60'000 && temp <= 200'000)
    cout << temp << " " << rad << '\n';


The rad variable used in the temperature equation in the for loop should use the rad value from the input.


If you remove the rad calculation on L36 then within the for loop nothing changes so you will obtain the same result for all values of k... Are you trying to use rad for 2 different purposes?
Last edited on
lastchance (6980)
Try this one.

It's based on the linear (for density) stellar model in https://www.astro.utu.fi/~cflynn/Stars/l3.html
and uses equation 29 (or 48) for central density and equation 34 (or 43) for radial distribution of temperature. (This appears to be what you are using).

NOTE: I have switched to mks (metre-kilogram-second) units, rather than the cgs (centimetre-gram-second) units in that file.

NOTE: I have used the mass of the star to get central density, but if you want to use the temperature at, e.g., the core then you can simply set T0 instead of computing it.

It's not the greatest of models: the absolute temperature is predicted to be 0 at the surface of the star and an ideal-gas approximation will be highly dubious at those plasma-inducing temperatures.

The test-case data below is for the sun. The range of radii for which the temperature would lie between 60 000 and 200 000 K is very small (corona only).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
 
#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;

int main()
{
   const double PI = 4.0 * atan( 1.0 );        // pi
   const double G = 6.6743e-11;                // universal gravitational constant in N m^2 / kg^2
   const double kb = 1.380649e-23;             // Boltzmann constant in J / K (EXACT, as of the 2019 redefinition of SI units)
   const double mp = 1.67262192e-27;           // proton mass in kg
   double mu = 0.6;                            // ratio of the average mass of subatomic particles in the star to that of a proton
  
   double R;                                   // star radius in m
   double M;                                   // star mass in kg
   int n;                                      // number of points to be plotted

   cout << "Enter the radius of the star in m: ";   cin >> R;
   cout << "Enter the mass of the star in kg:  ";   cin >> M;
   cout << "Enter number of radii (inc. ends): ";   cin >> n;

   double dr = R / ( n - 1 );
   double rho_centre = 3 * M / ( PI * R * R * R );    // central density (linear variation - see eqn 29)
   double T0 = ( 5.0 * PI * G * mu * mp ) / ( 36.0 * kb ) * rho_centre * R * R;
                                                      // core temperature (equals the constant in equation 34)

   #define FMT << setw( 15 ) << setprecision( 2 ) << fixed <<
   cout << "\nCore temperature (K): " FMT T0 << "\n\n";
   cout FMT "r/R" FMT "T(K)" << '\n';
   for ( int i = 1; i <= n; i++)
   { 
      double r = ( i - 1 ) * dr;                                               // radius (m)
      double rR = r / R;                                                       // relative radius
      double T = T0 * ( 1.0 + rR - 3.8 * rR * rR + 1.8  * rR * rR * rR );      // absolute temperature (K)
      cout FMT rR FMT T << '\n';
   } 
}


Enter the radius of the star in m: 6.963e8
Enter the mass of the star in kg:  1.989e30
Enter number of radii (inc. ends): 11

Core temperature (K):      5774290.40

            r/R           T(K)
           0.00     5774290.40
           0.10     6142690.13
           0.20     6134606.12
           0.30     5812400.72
           0.40     5238436.25
           0.50     4475075.06
           0.60     3584679.48
           0.70     2629611.85
           0.80     1672234.50
           0.90      774909.77
           1.00           0.00
Last edited on
Topic archived. No new replies allowed.