 |
 |
| View previous topic :: View next topic |
| Author |
Message |
Jerry I
Joined: 14 Sep 2003
Posts: 101
Location: Toronto, Ontario, Canada
|
| Trig functions math not accurate values obtained. |
 Posted: Wed Mar 26, 2025 1:14 pm |
 |
|
Using a PIC18F27K40 28pin Dip
Compiler Version 5.112. PCH
The following sample code will do a calculation of the sun position and angle for the provided date & time. I would not expect the precision to be perfect but is there any suggestions you may provide to make output values better. I have attached a complete program that will compile and run. I will also attach the original source code which also compiles and runs on a online C-compiler to see the accurate results.
[url]
https://www.programiz.com/c-programming/online-compiler/
[/url]
The original code used mostly the double type. I changed the type to float. There are fprintf output statements though out the code showing values. I was not expecting the %f printf usage function printing invalid calculations. This is shown in the code. The variable JD-1, JD-2.
This is the Original Code I converted to compile on the Pic.
Copy and paste into the online C-Compiler and run to see results.
| Code: |
#include <stdio.h>
#include <math.h>
#define PI 3.14159265358979323846
#define deg2rad(x) ((x) * PI / 180.0)
#define rad2deg(x) ((x) * 180.0 / PI)
double compute_jd(int year, int month, int day, int hour, int minute, double tz) {
// Adjust for time zone to get UTC
double utc_hour = hour - tz;
double utc_minute = minute;
printf("\nutc_hr = %lf utc_min %lf\n", utc_hour, utc_minute);
// Handle day overflow
int day_adjust = 0;
if (utc_hour < 0) { utc_hour += 24; day_adjust = -1; }
if (utc_hour >= 24) { utc_hour -= 24; day_adjust = 1; }
// Simplified Julian Date calculation
int a = (14 - month) / 12;
int y = year + 4800 - a;
int m = month + 12*a - 3;
printf("a = %d y = %d m = %d\n", a,y,m);
double jd = day + (153*m + 2)/5 + 365*y + y/4 - y/100 + y/400 - 32045;
printf("JD-1 = %lf\n", jd);
jd += (utc_hour + utc_minute/60.0) / 24.0 - 0.5;
printf("JD-2 = %lf\n\n", jd);
return jd;
}
void sun_position(double jd, double lat, double lon, double *azimuth, double *altitude) {
double T = (jd - 2451545.0) / 36525.0;
// Sun's coordinates
double M = fmod(357.52911 + 35999.05029*T, 360.0);
double M_rad = deg2rad(M);
double C = deg2rad(1.914602*sin(M_rad) + 0.019993*sin(2*M_rad) + 0.000289*sin(3*M_rad));
double L0 = deg2rad(fmod(280.46646 + 36000.76983*T, 360.0));
double L = L0 + C;
printf( "T = %f M = %f M_rad = %f C = %f L0 = %f L = %f\n\n",T,M,M_rad,C,L0,L);
//a - trunc(a/b) * b
// Ecliptic coordinates
double epsilon = deg2rad(23.43929111 - 0.0130042*T);
double RA = atan2(cos(epsilon)*sin(L), cos(L));
double delta = asin(sin(epsilon)*sin(L));
// Sidereal time
double GMST = fmod(280.46061837 + 360.98564736629*(jd - 2451545.0), 360.0);
double H = deg2rad(fmod(GMST + lon - rad2deg(RA) + 180.0, 360.0) - 180.0);
printf("RA->%f\r\n",rad2deg(RA));
printf("epsilon = %2.5f RA = %2.5f delta = %2.5f GMST = %2.5f H = %2.5f \n\n",epsilon,RA,delta,GMST,H);
// Convert to horizontal coordinates
double lat_rad = deg2rad(lat);
double sin_alt = sin(delta)*sin(lat_rad) + cos(delta)*cos(lat_rad)*cos(H);
*altitude = rad2deg(asin(sin_alt));
double cos_azm = (sin(delta) - sin(lat_rad)*sin_alt) / (cos(lat_rad)*cos(asin(sin_alt)));
*azimuth = rad2deg(acos(cos_azm));
if (sin(H) > 0) *azimuth = 360.0 - *azimuth;
}
double compute_tilt(double altitude, double azimuth) {
if (altitude <= 0) return 0.0;
double gamma = azimuth - 180.0; // Convert to south-based
double tilt = rad2deg(atan(sin(deg2rad(gamma)) / tan(deg2rad(altitude))));
return tilt + 0.0; // Add fixed 20-degree tilt
}
int main() {
int year, month, day, hour, minute;
double tz, lat, lon;
// printf("Enter year: "); scanf("%d", &year);
printf("Enter month: "); scanf("%d", &month);
printf("Enter day: "); scanf("%d", &day);
printf("Enter hour: "); scanf("%d", &hour);
printf("Enter minute: "); scanf("%d", &minute);
// printf("Enter time zone offset: "); scanf("%lf", &tz);
// printf("Enter latitude: "); scanf("%lf", &lat);
// printf("Enter longitude: "); scanf("%lf", &lon);
year = 2025;
lat = 43.6532;
lon = -79.3832;
tz = -5;
double jd = compute_jd(year, month, day, hour, minute, tz);
double azimuth, altitude;
sun_position(jd, lat, lon, &azimuth, &altitude);
double tilt = compute_tilt(altitude, azimuth);
printf("Calculating Sun Position for time %02d/%02d/%d %02d:%02d:00\n", month, day, year, hour, minute);
printf("GEO %lf, %lf\n", lat, lon);
printf("\nSun Position:\n");
printf(" Azimuth: %.2f° Altitude: %.2f°\n", azimuth, altitude);
printf(" East/West Optimal Tilt Angle from Horizontal: %.2f°\n\n", tilt);
return 0;
}
|
Online compiler output
Enter month: 6
Enter day: 21
Enter hour: 10
Enter minute: 45
utc_hr = 15.000000 utc_min 45.000000
a = 0 y = 6825 m = 3
JD-1 = 2460848.000000
JD-2 = 2460848.156250
T = 0.254707 M = 166.722531 M_rad = 2.909857 C = 0.007522 L0 = 1.572504 L = 1.580026
RA->90.576372
epsilon = 0.40903 RA = 1.58086 delta = 0.40902 GMST = 146.34207 H = -0.41220
Calculating Sun Position for time 06/21/2025 10:45:00
GEO 43.653200, -79.383200
Sun Position:
Azimuth: 128.52° Altitude: 61.98°
East/West Optimal Tilt Angle from Horizontal: -22.61°
Code output sample from Pic.
Enter month: 6
Enter day: 21
Enter hour: 10
Enter minute: 45
utc_hr = 15.00 utc_min 45.00
a = 0 y = 6825 m = 3
JD-1 %2.5f value = 12716.76939 -> ??? %f value = 2460848.00
JD-1 %2.5f value = 12716.76939 -> ??? %f value = 2460848.16
T = 0.25470 M = 166.81641 M_rad = 2.91149 C = 0.00746 L0 = 1.57413 L = 1.58160
RA->90.67506
epsilon = 0.40903 RA = 1.58257 delta = 0.40900 GMST = 180.25 H = 0.17
Calculating Sun Position for time 06/21/2025 10:45:00
GEO (43.653, -79.383)
Sun Position:
Azimuth: 205.81 Altitude: 68.11
East/West Optimal Tilt Angle from Horizontal: 29.92
Pic running test code below
| Code: |
#define PIC_18F27K40
#include <18F27K40.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#zero_ram
#zero_local_ram
#define BAUD1_SPEED 9600
#define BAUD2_SPEED 19200
#case
#fuses NOWDT, NOPPS1WAY
#fuses MCLR, PUT, BROWNOUT, BORV27, NOCKS, EBTR
// PIC18F27K40 PIN SELECTION
#PIN_SELECT U1RX=PIN_C7
#PIN_SELECT U1TX=PIN_C6
#define I2C_SDA PIN_A7
#define I2C_SCL PIN_A6
#use delay(INTERNAL, CLOCK=32MHz) /////////// Check Crystal Value ////////////
#use i2c(Master, Slow=100000, sda=I2C_SDA, scl=I2C_SCL, force_sw)
#use rs232(baud=9600, UART1, bits=8, parity=N, ERRORS, stream=HOSTPC)
#use FAST_IO(C)
#PRIORITY rda
#define PI 3.14159265358979323846
#define deg2rad(x) ((x) * PI / 180.0)
#define rad2deg(x) ((x) * 180.0 / PI)
float compute_jd(int16 year, int month, int day, int hour, int minute, signed int tz) {
// Adjust for time zone to get UTC
float utc_hour = hour - tz;
float utc_minute = minute;
fprintf(HOSTPC, "\r\nutc_hr = %f utc_min %f\r\n", utc_hour, utc_minute);
// Handle day overflow
int day_adjust = 0;
if (utc_hour < 0) { utc_hour += 24; day_adjust = -1; }
if (utc_hour >= 24) { utc_hour -= 24; day_adjust = 1; }
// Simplified Julian Date calculation
signed int32 a = (14 - month) / 12;
signed int32 y = year + 4800 - a;
signed int32 m = month + 12*a - 3;
fprintf(HOSTPC,"a = %Ld y = %Ld m = %Ld\r\n", a,y,m);
float jd = day + (153*m + 2)/5 + 365*y + y/4 - y/100 + y/400 - 32045;
fprintf(HOSTPC,"JD-1 %%2.5f value = %2.5f -> ??? %%f value = %f\r\n", jd, jd);
// ***** ERROR ABOVE AND BELOW for jd value DISPLAYED %2.5F *********
// The %2.5f value is way off when displayed.
// It looks like it gives the priority to the the # of decimal points
// I think that only 10 digits are allowed to be displayed ??
jd += (utc_hour + utc_minute/60.0) / 24.0 - 0.5;
fprintf(HOSTPC,"JD-1 %%2.5f value = %2.5f -> ??? %%f value = %f\r\n\r\n", jd, jd);
//fprintf(HOSTPC,"JD-2 = %2.5f -> ??? %f\r\n\r\n", jd, jd);
return jd;
}
void sun_position(float jd, float lat, float lon, float *azimuth, float *altitude) {
float T = (jd - 2451545.0) / 36525.0;
// Sun's coordinates
float M = (float)fmod(357.52911 + 35999.05029*T, 360.0);
float M_rad = deg2rad(M);
float C = (float)deg2rad(1.914602*sin(M_rad) + 0.019993*sin(2*M_rad) + 0.000289*sin(3*M_rad));
float L0 = (float)deg2rad(fmod(280.46646 + 36000.76983*T, 360.0));
float L = L0 + C;
//fprintf(HOSTPC, "T = %2.5f M = %2.5f M_rad = %2.5f C = %2.5f L0 = %f L = %f\r\n\r\n",T,M,M_rad,C,L0,L);
fprintf(HOSTPC, "T = %2.5f M = %2.5f M_rad = %2.5f C = %2.5f L0 = %2.5f L = %2.5f\r\n\r\n",T,M,M_rad,C,L0,L);
/*
// Ecliptic coordinates
float epsilon = (float)deg2rad(23.43929111 - 0.0130042*T);
float RA = atan2(cos(epsilon)*sin(L), cos(L));
float RA = atan2(cos(epsilon)*sin(L), cos(L));
float delta = asin(sin(epsilon)*sin(L));
*/
// Ecliptic coordinates
float epsilon = deg2rad(23.43929111 - 0.0130042*T);
float RA = atan2(cos(epsilon)*sin(L), cos(L));
float delta = asin(sin(epsilon)*sin(L));
//printf("epsilon = %2.5f RA = %2.5f delta = %2.5f\r\n",epsilon,RA,delta);
float GMST = fmod(280.46061837 + 360.98564736629*(jd - 2451545.0), 360.0);
float H = (float)deg2rad(fmod(GMST + lon - (float)rad2deg(RA) + 180.0, 360.0) - 180.0);
fprintf(HOSTPC, "RA->%2.5f\r\n\r",rad2deg(RA));
fprintf(HOSTPC, "epsilon = %2.5f RA = %2.5f delta = %2.5f GMST = %f H = %f \r\n\r\n",epsilon,RA,delta,GMST,H);
// Convert to horizontal coordinates
float lat_rad = deg2rad(lat);
float sin_alt = sin(delta)*sin(lat_rad) + cos(delta)*cos(lat_rad)*cos(H);
*altitude = rad2deg(asin(sin_alt));
float cos_azm = (sin(delta) - sin(lat_rad)*sin_alt) / (cos(lat_rad)*cos(asin(sin_alt)));
*azimuth = rad2deg(acos(cos_azm));
if (sin(H) > 0) *azimuth = 360.0 - *azimuth;
}
float compute_tilt(float altitude, float azimuth) {
if (altitude <= 0) return 0.0;
float gamma = azimuth - 180.0; // Convert to south-based
float tilt = rad2deg(atan(sin(deg2rad(gamma)) / tan(deg2rad(altitude))));
return tilt + 20.0; // Add fixed 20-degree tilt
}
int main() {
int month, day, hour, minute;
int16 year;
float tz, lat, lon;
// printf("Enter year: "); scanf("%d", &year);
fprintf(HOSTPC, "\r\n\r\nEnter month: "); fscanf(HOSTPC, "%d", &month); fprintf(HOSTPC, "%d\r\n", month);
fprintf(HOSTPC, "Enter day: "); fscanf(HOSTPC, "%d", &day); fprintf(HOSTPC, "%d\r\n", day);
fprintf(HOSTPC, "Enter hour: "); fscanf(HOSTPC, "%d", &hour); fprintf(HOSTPC, "%d\r\n", hour);
fprintf(HOSTPC, "Enter minute: "); fscanf(HOSTPC, "%d", &minute); fprintf(HOSTPC, "%d\r\n", minute);
// printf("Enter time zone offset: "); scanf("%lf", &tz);
// printf("Enter latitude: "); scanf("%lf", &lat);
// printf("Enter longitude: "); scanf("%lf", &lon);
year = 2025;
lat = 43.6532;
lon = -79.3832;
tz = -5;
float jd = compute_jd(year, month, day, hour, minute, tz);
float azimuth, altitude;
sun_position(jd, lat, lon, &azimuth, &altitude);
float tilt = compute_tilt(altitude, azimuth);
fprintf(HOSTPC, "Calculating Sun Position for time %02d/%02d/%ld %02d:%02d:00\r\n", month, day, year, hour, minute);
fprintf(HOSTPC, "GEO (%2.3f, %2.3f)\r\n", lat, lon);
fprintf(HOSTPC, "\r\nSun Position:\r\n");
fprintf(HOSTPC, " Azimuth: %.2f Altitude: %.2f\r\n", azimuth, altitude);
fprintf(HOSTPC, " East/West Optimal Tilt Angle from Horizontal: %.2f\r\n", tilt);
//°
reset_cpu();
return 0;
}
|
I am sorry for the long post, but it is easier for someone to try the code to see the calculation errors. From the output values given the GMST & H variables are the most incorrect.
Thank you again for any help you can provide. |
|
 |
temtronic
Joined: 01 Jul 2010
Posts: 9442
Location: Greensville,Ontario
|
|
| Posted: Wed Mar 26, 2025 2:03 pm |
|
|
Hmm, I'll try to find my code for a 'sunrise/set chicken coop project from a few years ago. I remember 'ansi' was part of the solution........
The 64K question is... WHICH PC was a using back then ?????
surprise I found them but they figure out sunrise and sunset, don't know if it's what you need ? |
|
|
Jerry I
Joined: 14 Sep 2003
Posts: 101
Location: Toronto, Ontario, Canada
|
|
| Posted: Wed Mar 26, 2025 2:24 pm |
|
|
Thank you for info.
This code calculates the tilt angle of a 1x solar tracker that tilts East <-> West on a North South axis. The face of the solar panels can also be tilted down to the south to attract more sunlight. |
|
|
Jerson
Joined: 31 Jul 2009
Posts: 127
Location: Bombay, India
|
|
| Posted: Wed Mar 26, 2025 8:50 pm |
|
|
I'm not sure if this hint will help you in any way. AFAIK - CCS-C does not have a double precision type or perhaps it is implemented as single internally - please check on that.
Many years ago, I did a similar calculation on a 8051 processor design. I do recall vividly that the compiler needs to have 'double precision' floating point support for the astronomical calculations to be near precise. 'single precision' floats do not cut it.
I did a side-by-side 'difference' of your 2 codes and they look alike. Therefore, the reason cited above could be why you have this problem you state. |
|
|
Ttelmah
Joined: 11 Mar 2010
Posts: 19755
|
|
| Posted: Thu Mar 27, 2025 1:54 am |
|
|
Double is only available with the PIC24/30/33 chips.
Normally the errors between double and single are not enough to matter.
I did calculations of planet position using single precision, and was within a
tiny fraction of a degree of the real value. The errors being shown here
are much larger than the maths errors should give, suggests there is
something wrong in the code somewhere. I would try splitting some
of the long maths parts up. CCS C has problems if too much arithmetic
is done in a single line. Can result in internal overflows that it doesn't
handle. There are several lines in the code that might well trigger this.
So do the maths a part at a time, and check the result after each piece.
The trig functions generally can be improved a little if you are prepared to
take more time. Also I did code a library for myself in the past using
BCD arithmetic. Bulky, and very slow, but you could achieve much better
precision.
Choosing algorithms to be less affected by maths errors is an important
exercise sometimes.
Try an experiment. Run the algorithm using float32 maths online. If the
result differs from the PIC version, then step through the maths section
by section, till you find what is going wrong.
Wonder if somebody has the time to try on a current compiler?. There
was a version with a maths error around your version. Might be the
problem. Jerson says he has dome this?.
Remembering back, it was the atan2 function that had an issue. There
were posts about this here at the time.
The maths issue, was corruption of various registers when arithmetic
functions were combined with using printf:
[url]
https://www.ccsinfo.com/forum/viewtopic.php?t=51298&highlight=atan2
[/url]
It is a really horrible insidious problem giving different values according to
what you are doing. |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ 
Search
Register
Profile 
Log in to check your private messages 
Log in
