Problem with exp() function in C?

[vboro]

Greetings!

I have a C code in which I count the descreasing of a satellite. I translated it to C language from the original (QBasic, see link below) but the exp() function doesn't seem to work for some reason. Any clues how I could fix it?

Thanks for the help in advance!!

My translated C code:

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

int main(){

  double ms, sa, h, rf, ga;

  printf("Mass of satellite (kg): ");
  scanf("%lf", &ms);

  printf("A of satellite (m^2): ");
  scanf("%lf", &sa);

  printf("Distance of satellite from Earth (km): ");
  scanf("%lf", &h);

  printf("Solar Radio Flux: ");
  scanf("%lf", &rf);

  printf("Geomagnetic A index: ");
  scanf("%lf", &ga);

  double Me = 5.98e+24, G = 6.67e-11, pi = 3.1416;
  int T = 0, dT = 1, H1 = 10, H2 = h;
  long Re = 6378000;
  long D9 = dT * 3600 * 24, R = Re + h * 1000;
  double P = 2 * pi * pow(R * R * R / Me / G, 2);

  printf("Time\t\tHeight\t\tPeriod\t\tMean motion\t\tDecay");
  
  do{
    int SH = (900 + 2.5 * (rf - 70 ) + 1.5 * ga) / ( 27 - (3/250)*(h-200));
    double DN = (6e-10)*exp(-((h-175)/SH)); //atmospheric density
    int dP = 3*pi*sa/ms*R*DN*D9; //decrement in orbital period
    if(h <= H2){ //test for print
      printf("h is lower than H2!\n");
      int Pm = P / 60;
      double MM = 1440 / Pm; 
      double nMM = 1440 / ((P - dP)) / 60; //print units
      int Decay = dP / dT / P * MM ; //rev/day/day
      //LPRINT USING f$ ; T ; H ; P / 60 ; MM ; Decay ;
      P = P/60;
      printf("%d \t %lf \t %lf \t %lf \t %d \n", T, h, P, MM, Decay);
//do print
      H2 = H2 - H1; //decrement print height
    }
    //end of print routine
    if( h < 180 ){ 
      break; 
    }
    
      P = P - dP; 
      T = T + dT; //compute new values
      R = pow((G * Me * P * P / 4 / pi / pi ), (1/3)); //new orbital radius
      h = (R - Re)/1000; //new altitude (semimajor axis)

    }while(h >= 180);

    printf("Re-entry after %d days", T);
    
  return 0;
}

The original code (page 4): https://www.sws.bom.gov.au/Category/Educational/Space%20Weather/Space%20Weather%20Effects/SatelliteOrbitalDecayCalculations.pdf

[lastchance]

You have integer division problems all over the place (including inside your exp function ... but also elsewhere).

3 /250 -> 0 with integer division (line 33)
(h-175)/SH is going to do integer division and get it wrong
1/3 -> 0 with integer division.

[Peter87]

lastchance wrote:

(h-175)/SH is going to do integer division and get it wrong

This particular expression is not an issue since h is a double.

[vboro]

Okay, I changed the integers to doubles :) But the exp() function is still a problem... any ideas what can the problem be?

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

int main(){

  double ms, sa, h, rf, ga;

  printf("Mass of satellite (kg): ");
  scanf("%lf", &ms);

  printf("A of satellite (m^2): ");
  scanf("%lf", &sa);

  printf("Distance of satellite from Earth (km): ");
  scanf("%lf", &h);

  printf("Solar Radio Flux: ");
  scanf("%lf", &rf);

  printf("Geomagnetic A index: ");
  scanf("%lf", &ga);

  double Me = 5.98e+24, G = 6.67e-11, pi = 3.1416;
  double T = 0, dT = 1, H1 = 10, H2 = h;
  long Re = 6378000;
  long D9 = dT * 3600 * 24, R = Re + h * 1000;
  double P = 2 * pi * pow(R * R * R / Me / G, 2);

  do{
    double SH = (900 + 2.5 * (rf - 70 ) + 1.5 * ga) / ( 27 - (3/250)*(h-200));
    double DN = (6e-10)*exp(-((h-175)/SH)); //atmospheric density
    double dP = 3*pi*sa/ms*R*DN*D9; //decrement in orbital period
    if(h <= H2){ //test for prdouble
      //prdoublef("", );
      double Pm = P / 60;
      double MM = 1440 / Pm; 
      double nMM = 1440 / ((P - dP)) / 60; //prdouble units
      double Decay = dP / dT / P * MM ; //rev/day/day
      printf("%lf \t %lf \t %lf \t %lf \t %lf \n", T, h, P, MM, Decay);
      H2 = H2 - H1; //decrement prdouble height
    }
    //end of prdouble routine
    if( h < 180 ){ 
      break; 
    }
    
      P = P - dP; 
      T = T + dT; //compute new values
      R = pow((G * Me * P * P / 4 / pi / pi ), (1/3)); //new orbital radius
      h = (R - Re)/1000; //new altitude (semimajor axis)


    }while(h < 180);

  printf("Re-entry after %lf days", T);
    
  return 0;
}

[lastchance]

@vboro, you ignored my other points. Both 3/250 and 1/3 still give 0 by integer division in C++.

Try this:

#include <iostream>
#include <cmath>
using namespace std;

int main()
{
   const double Me = 5.98e+24;
   const double Re = 6378000;
   const double G = 6.67e-11;
   const double pi = 4.0 * atan( 1.0 );

   double ms;
   cout << "Mass of satellite (kg): ";   cin >> ms;
  
   double sa;
   cout << "A of satellite (m^2): ";   cin >> sa;  
  
   double h;
   cout << "Distance of satellite from Earth (km): ";   cin >> h;     
  
   double rf;
   cout << "Solar Radio Flux: ";   cin >> rf;
  
   double ga;
   cout << "Geomagnetic A index: ";   cin >> ga;


   double T = 0.0, dT = 0.1, H1 = 10.0, H2 = h;
   double D9 = dT * 3600 * 24;
   double R = Re + h * 1000;
   double P = 2 * pi * sqrt( R * R * R / Me / G );


   cout << "Time\t\tHeight\t\tPeriod\t\tMean motion\t\tDecay\n";
   while( true )
   {
      double SH = ( 900 + 2.5 * ( rf - 70 ) + 1.5 * ga ) / ( 27 - 0.012 * ( h - 200 ) );
      double DN = 6.0e-10 * exp( -( h - 175 ) / SH );      // atmospheric density
      double dP = 3 * pi * ( sa / ms ) * R * DN * D9;      // decrement in orbital period

      if ( h <= H2 )
      {
          double Pm = P / 60;
          double MM = 1440 / Pm;
          double Decay = dP / dT / P * MM ; //rev/day/day
          cout << T << "\t\t" << h << "\t\t" << P << "\t\t" << MM << "\t\t" << Decay << '\n';
          H2 -= H1; //decrement print height
      }
      if( h < 180 ) break;
   
      P -= dP;
      T += dT;
      R = pow( G * Me * P * P / 4 / pi / pi, 1.0 / 3.0 );  // new orbital radius
      h = ( R - Re ) / 1000;                               // new altitude (semimajor axis)

   }

   cout << "Re-entry after " << T << " days";
}

Mass of satellite (kg): 100
A of satellite (m^2): 1
Distance of satellite from Earth (km): 300
Solar Radio Flux: 70
Geomagnetic A index: 0
Time		Height		Period		Mean motion		Decay
0		300		5429.2		15.9139		0.00265707
11.9		289.905		5416.9		15.9501		0.00350492
20.8		279.999		5404.83		15.9857		0.00461164
27.7		269.898		5392.53		16.0222		0.00611716
32.9		259.838		5380.3		16.0586		0.00812687
36.8		249.859		5368.17		16.0949		0.010801
39.8		239.665		5355.79		16.1321		0.0144828
42		229.716		5343.72		16.1685		0.0193347
43.7		219.458		5331.28		16.2062		0.0261175
44.9		209.814		5319.59		16.2418		0.0347396
45.9		199.022		5306.53		16.2818		0.0479446
46.6		188.762		5294.11		16.32		0.065316
47.1		178.94		5282.24		16.3567		0.088044
Re-entry after 47.1 days 

[vboro]

Ooo so those numbers needed a change, as well ^^' Sorry, I missed that part, I was trying to act as you said. But many thanks for the final corrections, I can't be more thankful!!

[George P]

C++20 added several mathematical constants to the language, including pi.

#include <iostream>
#include <numbers>  // https://en.cppreference.com/w/cpp/numeric/constants

int main()
{
   std::cout << std::numbers::pi << '\n';
}

[TheIdeasMan]

As an aside, in C programming with scanf, one should check the value returned, to see if it worked or not. This is mainly due to not knowing whether the input is in the correct format or not. The other thing is assumptions about consumption of white space can cause errors.

https://en.cppreference.com/w/c/io/fscanf

Note that all the versions return an int

One should use the secure versions (with the _s suffix) to avoid problems like buffer overflow attacks.

If one wanted to be thoroughly pedantic, one could check the return value from printf as well.

In C++ std::printf is different from C's printf in that it is statically type safe. That is in the standard somewhere.

There are third party C++ formatted input/output library functions that are dynamically type and thread safe, with similar formatting options to printf/scanf. They are useful when there is a lot of format changes, one can avoid pages of normal C++ code, in favour of more dense printf like code. For example, try coding this with std::cout, all on 1 line of output, compare with printf style code:

All fields right justified

table boundary | ;
line number, 2 digits padded with 0;
field separator | ;
degrees, 3 digits, 0 padded;
degree symbol;
minutes, 2 digits, 0 padded;
minutes symbol ' ;
seconds, 2 digits, 0 padded;
seconds symbol " ;
field separator | ;
Distance width 9, 3 dp;
field separator | ;
Easting width 9, 3 dp;
field separator | ;
Easting adjustment, 2 digits with - sign if negative;
field separator | ;
Northing width 9, 3 dp;
field separator | ;
Northing adjustment, 2 digits with - sign if negative;
field separator | ;
EAST width 11, 3 dp;
field separator | ;
NORTH width 11, 3 dp;
Comment/ Description - as many chars available (left over) to fit when printed A4 landscape
table boundary | ;

[George P]

TheIdeasMan wrote:

One should use the secure versions (with the _s suffix) to avoid problems like buffer overflow attacks.

Not too long ago I was reminded not all implementations provide the bounds checking functions that C11 included in the standard.

http://www.cplusplus.com/forum/beginner/282327/#msg1221737

Yes, I totally agree one should use the bounds checking functions when available.

From https://en.cppreference.com/w/c/io/fscanf:

As with all bounds-checked functions, scanf_s , fscanf_s, and sscanf_s are only guaranteed to be available if __STDC_LIB_EXT1__ is defined by the implementation and if the user defines __STDC_WANT_LIB_EXT1__ to the integer constant 1 before including stdio.h.

If using MSVC you get the bounds checking functions "for free" without any #define.

[Peter87]

In other words, don't use the _s versions if you're writing standard C++, or in C if you want it to work with all implementations.

[seeplus]

There are third party C++ formatted input/output library functions that are dynamically type and thread safe, with similar formatting options to printf/scanf

std::format et al in C++20.
https://en.cppreference.com/w/cpp/utility/format

[TheIdeasMan]

@seeplus

With std::format, I should have mentioned that, it looks easier. However it is not implemented yet in g++ or clang++, only in v19.29 MSVC.

[George P]

There is the fmt library if std::format isn't available. The API looks remarkably similar to what the C++ standard uses.

https://fmt.dev/latest/index.html

[TheIdeasMan]

George P wrote:

The API looks remarkably similar to what the C++ standard uses.

Yep, std::format was modeled on fmt :+)

[Ganado]

It became part of the borg.

[George P]

TheIdeasMan wrote:

Yep, std::format was modeled on fmt :+)

Of that I have no doubts. Wouldn't surprise me if the standards committee "stole" it lock, stock and barrel as they do with so many parts of the Boost libraries.

The fmt library has one part that is not part of std::format AFAIK, the fmt::print() function. Being able to draw basic graphics has its uses:

fmt::print(
  "┌{0:─^{2}}┐\n"
  "│{1: ^{2}}│\n"
  "└{0:─^{2}}┘\n", "", "Hello, world!", 20);

From https://fmt.dev/latest/syntax.html (Box drawing using Unicode fill)

[George P]

Coincidentally I was looking over Boost's format library and I noticed how very similar the formatting syntax was to C's formatting syntax. With a C++ twist, of course.

C++20 std::format/fmt library took a fork in the road route.

[Peter87]

I've heard people say fmt/std::format was inspired by Python. I don't use Python so I'm not sure but I think at least the {} syntax was inspired by Python.

[George P]

https://en.cppreference.com/w/cpp/utility/format/formatter#Standard_format_specification

For basic types and string types, the format specification is based on the format specification in Python

Here's a bit less technical explanation of the formatting for "us dummies" just now learning about std::format:
https://docs.python.org/3/library/string.html#formatspec