Precision calculation problems in arduino DUE

sagise · February 8, 2022, 2:45pm

I have an arduino DUE to make a program that requires high-precision calculations, but the results are not as expected by comparing the results that I calculated using Excel which has tested the results.
I include screenshots and sketches to make it easier for you to understand.

I hope someone can help with my project, thanks for those who want to think about this.

 #include <math.h>
 int l=-6,b=109,tgl=13,bln=6,thn=2020,j=1,mnt=0,dtk=0,zw=(b/15);
 
      int M,Y,A,B;
    
      double Tahun = thn+(bln-1)/12+tgl/365;
    
      double delta;
      double bujurRataM;
      double anomaliRataM;
      double koreksiM;
      double bujurEkliptikaM;
      double omegaM;
      double epsilonZero;
      double deltaEpsilon;
      double epsilon;
      double JDUT;
      double TuJD;
      double GST_UT;
      double GST_L;
      double LST;
      double bujurNampakM;
      double alphaM;
      double deltaM;
      double hourAngle;
      double azimuthMselatan;
      double azimuthM;
      double altitudM;
     
    void setup()
    {
      
      Serial.begin (115200);
      
     
      if (bln<3) 
      M = bln+1;
    else 
      M = bln;
    
    if (bln<3) 
      Y = thn-1;
    else 
      Y = thn;
    
    A=(Y/100);
    
    B = (2+int (A/4)-A);
    
    if (Tahun>2005)
      if(Tahun<=2050)
        delta = 62.92 + 0.32217*(Tahun-2000) + 0.005589*(Tahun-2000)*(Tahun-2000);
      else 
        delta = 0; 
    else
      delta = 0;
    
  double JDU = 1720994.5+int32_t(365.25*Y)+int32_t(30.60001*(M+1))+B+tgl+(j+mnt/60+dtk/3600)/24 - zw/24;
    double JDE = JDU+(delta/86400);
      double T = (JDE-2451545)/36525;
    
    
    JDUT = 1720994.5+int32_t(365.25*Y)+int32_t(30.60001*(M+1))+B+j;
      TuJD = (JDUT-2451545)/36525;  
        GST_UT = fmod(6.6973745583+2400.0513369072*TuJD+0.0000258622*TuJD*TuJD,24);
          GST_L = fmod(GST_UT+(j+mnt/60+dtk/3600)*1.00273790935,24);
            LST = fmod(GST_L+b/15,24);
      
    //Posisi Matahari
    bujurRataM = fmod(280.46646+36000.76983*T+0.0003032*T*T,360);
      anomaliRataM = fmod(357.5291092 + 35999.0502909*T - 0.0001536*T*T + T*T*T/24490000, 360);
        koreksiM = (1.914602 - 0.004817*bujurRataM - 0.000014*bujurRataM*bujurRataM)*sin(radians(anomaliRataM)) + (0.019993 - 0.000101*bujurRataM)*sin(2*radians(anomaliRataM)) + 0.000289*sin(3*radians(anomaliRataM));
          bujurEkliptikaM = bujurRataM + koreksiM;
    omegaM = fmod(125.04-1934.136*T,360);
    epsilonZero =  23 + 26/60 + 21.448/3600 - 46.815*T/3600 - 0.00059*T*T/3600 + 0.001813*T*T*T/3600;
      deltaEpsilon = 0.00256*cos(radians(omegaM));
        epsilon = epsilonZero + deltaEpsilon;
    bujurNampakM = bujurEkliptikaM-0.00569-0.00478*sin(radians(omegaM));
      alphaM = fmod(degrees(atan2(cos(radians(bujurNampakM)),cos(radians(epsilon))*sin(radians(bujurNampakM)))),360);
    deltaM = -degrees(asin(sin(radians(epsilon))*sin(radians(bujurNampakM))));  
        hourAngle = fmod(LST-alphaM,24)*15;
          azimuthMselatan = atan2(cos(radians(hourAngle))*sin(radians(l))-tan(radians(deltaM))*cos(radians(l)),sin(radians(hourAngle)));
            azimuthM = fmod(radians(azimuthMselatan)+180,360);
          altitudM = degrees(asin(sin(radians(l))*sin(radians(deltaM))+cos(radians(l))*cos(radians(deltaM))*cos(radians(hourAngle))));

          Serial.print("Y = " ); 
        Serial.println(Y);
        Serial.print("A = " );
        Serial.println(A);
        Serial.print("B = " );
        Serial.println(B);
        Serial.print("delta = " );
        Serial.println(delta);
        Serial.print("JDU = " );
        Serial.println(JDU);
        Serial.print("JDE = " );
        Serial.println(JDE);
        Serial.print("T = " );
        Serial.println(T);
        Serial.print("bujurRataM = " );
        Serial.println(bujurRataM);
        Serial.print("anomaliRataM = " );
        Serial.println(anomaliRataM);
        Serial.print("koreksiM = " );
        Serial.println(koreksiM);
        Serial.print("bujurEkliptikaM = " );
        Serial.println(bujurEkliptikaM);
        Serial.print("omegaM = " );
        Serial.println(omegaM);
        Serial.print("epsilonZero = " );
        Serial.println(epsilonZero);
        Serial.print("deltaEpsilon = " );
        Serial.println(deltaEpsilon);
        Serial.print("epsilon = " );
        Serial.println(epsilon);
        Serial.print("JDUT = " );
        Serial.println(JDUT);
        Serial.print("TuJD = " );
        Serial.println(TuJD);
        Serial.print("GST_UT = " );
        Serial.println(GST_UT);
        Serial.print("GST_L = " );
        Serial.println(GST_L);
        Serial.print("LST = " );
        Serial.println(LST);
        Serial.print("bujurNampakM = " );
        Serial.println(bujurNampakM);
        Serial.print("alphaM = " );
        Serial.println(alphaM);
        Serial.print("deltaM = " );
        Serial.println(deltaM);
        Serial.print("hourAngle = " );
        Serial.println(hourAngle);
        Serial.print("azimuthMselatan = " );
        Serial.println(azimuthMselatan);
        Serial.print("azimuthM = " );
        Serial.println(azimuthM);
        Serial.print("altitudM = " );
        Serial.println(altitudM);


         }
    
    void loop()
    {
    }

gfvalvo · February 8, 2022, 2:52pm

Try printing the value of 'Tahun' and see if it's what you expect.

alto777 · February 8, 2022, 2:53pm

Is this just a matter of printing more decimal places?

DUE doubles are 8 bytes, so shoukd be accurate if you’ve transcribed the formulae correctly. Try

        Serial.println(alphaM, 6);

HTH

a7

anon73444976 · February 8, 2022, 3:00pm

Please could
you get
used to
using the auto-
format
tool before
posting your

code?

david_2018 · February 8, 2022, 3:11pm

Can you post the excel formula for the Julian date? There are parts of that equation that need to be done as integer math (dropping the fractional part). The excel value ending in something other than 0.50 looks suspicious.

sagise · February 8, 2022, 3:14pm

this is not only a matter of printing, but about the precision of the stored value results which will then be used for the next calculation.

sagise · February 8, 2022, 3:29pm

The formula I use in Excel is of course the same as the one in sketch

jremington · February 8, 2022, 3:29pm

The formula I use in Excel is of course the same as the one in sketch

Only if you wrote the Excel code correctly, and you gave us no reason to believe that you did.

sagise · February 8, 2022, 3:29pm

sorry, I don't understand what you mean

anon73444976 · February 8, 2022, 3:30pm

Your code wanders over the page like it is drunk, and makes it hard to read.
Using the IDE's auto-format tool would help fix that.

david_2018 · February 8, 2022, 3:59pm

Still looks wrong, Julian date should end in 0.50 because the day starts at noon instead of midnight.

Delta also looks incorrect, using a calculator gives a result of 71.599, which is what I get from the sketch you previously posted a few weeks ago with the same formula you are currently using (which originally comes from a NASA paper, I believe).

< edit > Also, it is very easy to calculate the Julian date from the unit timestamp, much simpler than doing the normal calculation, and usable as long as the year is after 1970.

alto777 · February 8, 2022, 4:02pm

It looks like Arduino doubles on DUE are according to the

IEEE 754 binary64 standard

Excel also adheres to IEEE 754, but I do not know what number of bytes are used.

I think you can expect to get similar or at least acceptable calculations.

TBH squinting at your screen shotz is without my capability. It would be easier to see what you are complaining about if you posted text we can all cut and paste and work with in our habitual ways.

a7

sagise · February 8, 2022, 4:12pm

does that mean I need to crop the screenshot to make it look clearer?

sagise · February 8, 2022, 4:15pm

can you give me the julian day calculation you mean?

gfvalvo · February 8, 2022, 4:15pm

OMG, it's not that difficult!!!

anon73444976 · February 8, 2022, 4:15pm

The only screenshots are the spreadsheet and the results, but I'm talking about your Arduino sketch.

sagise · February 8, 2022, 4:18pm

[code]



int l = -6, b = 109, tgl = 13, bln = 6, thn = 2020, j = 1, mnt = 0, dtk = 0, zw = (b / 15);

int M, Y, A, B;

double Tahun = thn + (bln - 1) / 12 + tgl / 365;

double delta;
double bujurRataM;
double anomaliRataM;
double koreksiM;
double bujurEkliptikaM;
double omegaM;
double epsilonZero;
double deltaEpsilon;
double epsilon;
double JDUT;
double TuJD;
double GST_UT;
double GST_L;
double LST;
double bujurNampakM;
double alphaM;
double deltaM;
double hourAngle;
double azimuthMselatan;
double azimuthM;
double altitudM;




void setup()
{

  Serial.begin (115200);


  if (bln < 3)
    M = bln + 1;
  else
    M = bln;

  if (bln < 3)
    Y = thn - 1;
  else
    Y = thn;

  A = (Y / 100);

  B = (2 + int (A / 4) - A);

  if (Tahun > 2005)
    if (Tahun <= 2050)
      delta = 62.92 + 0.32217 * (Tahun - 2000) + 0.005589 * (Tahun - 2000) * (Tahun - 2000);
    else
      delta = 0;
  else
    delta = 0;

  double JDU = 1720994.5 + int32_t(365.25 * Y) + int32_t(30.60001 * (M + 1)) + B + tgl + (j + mnt / 60 + dtk / 3600) / 24 - zw / 24;
  double JDE = JDU + (delta / 86400);
  double T = (JDE - 2451545) / 36525;


  JDUT = 1720994.5 + int32_t(365.25 * Y) + int32_t(30.60001 * (M + 1)) + B + j;
  TuJD = (JDUT - 2451545) / 36525;
  GST_UT = fmod(6.6973745583 + 2400.0513369072 * TuJD + 0.0000258622 * TuJD * TuJD, 24);
  GST_L = fmod(GST_UT + (j + mnt / 60 + dtk / 3600) * 1.00273790935, 24);
  LST = fmod(GST_L + b / 15, 24);

  //Posisi Matahari
  bujurRataM = fmod(280.46646 + 36000.76983 * T + 0.0003032 * T * T, 360);
  anomaliRataM = fmod(357.5291092 + 35999.0502909 * T - 0.0001536 * T * T + T * T * T / 24490000, 360);
  koreksiM = (1.914602 - 0.004817 * bujurRataM - 0.000014 * bujurRataM * bujurRataM) * sin(radians(anomaliRataM)) + (0.019993 - 0.000101 * bujurRataM) * sin(2 * radians(anomaliRataM)) + 0.000289 * sin(3 * radians(anomaliRataM));
  bujurEkliptikaM = bujurRataM + koreksiM;
  omegaM = fmod(125.04 - 1934.136 * T, 360);
  epsilonZero =  23 + 26 / 60 + 21.448 / 3600 - 46.815 * T / 3600 - 0.00059 * T * T / 3600 + 0.001813 * T * T * T / 3600;
  deltaEpsilon = 0.00256 * cos(radians(omegaM));
  epsilon = epsilonZero + deltaEpsilon;
  bujurNampakM = bujurEkliptikaM - 0.00569 - 0.00478 * sin(radians(omegaM));
  alphaM = fmod(degrees(atan2(cos(radians(bujurNampakM)), cos(radians(epsilon)) * sin(radians(bujurNampakM)))), 360);
  deltaM = -degrees(asin(sin(radians(epsilon)) * sin(radians(bujurNampakM))));
  hourAngle = fmod(LST - alphaM, 24) * 15;
  azimuthMselatan = atan2(cos(radians(hourAngle)) * sin(radians(l)) - tan(radians(deltaM)) * cos(radians(l)), sin(radians(hourAngle)));
  azimuthM = fmod(radians(azimuthMselatan) + 180, 360);
  altitudM = degrees(asin(sin(radians(l)) * sin(radians(deltaM)) + cos(radians(l)) * cos(radians(deltaM)) * cos(radians(hourAngle))));

  Serial.print(Tahun);
  Serial.print("Y = " );
  Serial.println(Y);
  Serial.print("A = " );
  Serial.println(A);
  Serial.print("B = " );
  Serial.println(B);
  Serial.print("delta = " );
  Serial.println(delta);
  Serial.print("JDU = " );
  Serial.println(JDU);
  Serial.print("JDE = " );
  Serial.println(JDE);
  Serial.print("T = " );
  Serial.println(T);
  Serial.print("bujurRataM = " );
  Serial.println(bujurRataM);
  Serial.print("anomaliRataM = " );
  Serial.println(anomaliRataM);
  Serial.print("koreksiM = " );
  Serial.println(koreksiM);
  Serial.print("bujurEkliptikaM = " );
  Serial.println(bujurEkliptikaM);
  Serial.print("omegaM = " );
  Serial.println(omegaM);
  Serial.print("epsilonZero = " );
  Serial.println(epsilonZero);
  Serial.print("deltaEpsilon = " );
  Serial.println(deltaEpsilon);
  Serial.print("epsilon = " );
  Serial.println(epsilon);
  Serial.print("JDUT = " );
  Serial.println(JDUT);
  Serial.print("TuJD = " );
  Serial.println(TuJD);
  Serial.print("GST_UT = " );
  Serial.println(GST_UT);
  Serial.print("GST_L = " );
  Serial.println(GST_L);
  Serial.print("LST = " );
  Serial.println(LST);
  Serial.print("bujurNampakM = " );
  Serial.println(bujurNampakM);
  Serial.print("alphaM = " );
  Serial.println(alphaM);
  Serial.print("deltaM = " );
  Serial.println(deltaM);
  Serial.print("hourAngle = " );
  Serial.println(hourAngle);
  Serial.print("azimuthMselatan = " );
  Serial.println(azimuthMselatan);
  Serial.print("azimuthM = " );
  Serial.println(azimuthM);
  Serial.print("altitudM = " );
  Serial.println(altitudM);


}

void loop()
{
}
[/code]

anon73444976 · February 8, 2022, 4:21pm

These are all integers.
See reply #2

gfvalvo · February 8, 2022, 4:23pm

I give advise once. If it's ignored, I don't repeat it.

sagise · February 8, 2022, 4:27pm

oke, then

Topic		Replies	Views
Decimal Precision Syntax & Programs	21	51513	May 6, 2021
[Résolu] Fonctions trigonométriques et variables à virgule flottante Français	36	11185	May 6, 2021
How does Ardupilot do it?? GPS accuracy issues !!! Distance Calculations Sensors	17	14336	May 6, 2021
Solved // Datum von gestern ohne Library berechnen Deutsch	22	10123	May 6, 2021
Mond berechnen… Deutsch	20	5588	May 6, 2021

Precision calculation problems in arduino DUE

Related topics