 |
 |
| View previous topic :: View next topic |
| Author |
Message |
Jerry I
Joined: 14 Sep 2003
Posts: 102
Location: Toronto, Ontario, Canada
|
| Trig functions math not accurate values obtained. UPDATE |
 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.
Last edited by Jerry I on Sat Mar 29, 2025 5:42 pm; edited 1 time in total |
|
 |
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: 102
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: 19761
|
|
| 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. |
|
|
Jerry I
Joined: 14 Sep 2003
Posts: 102
Location: Toronto, Ontario, Canada
|
|
| Posted: Sat Mar 29, 2025 5:33 pm |
|
|
CCS PCD C Compiler, Version 5.112, 59247 29-Mar-25 18:51
PIC24HJ256GP210
Thanks for any help. This is a real weird one. I have been pulling out my hair and I had very little to start with.
Ttelmah I tried what you asked, and to my surprise the online C-compiler also gave errors when I switched from double to floats. I also tried executing partial calculations but still does not work. Using floats in the declaration does not work for the Pic or the online C-compiler.
Ok now I am going nuts. I had a Pic24 development board so I compiled the same attached program using all variables that were float -> double. Everything compiles fine. But ????? I get no printf statements or any statements to not appear to execute in the code in the main() after the calls for the calculations I need. The delay_ms(5000) does not even work. The last statement is reset_cpu() and it does that.
Now if I change the all variables back to float (double -> float) the program compiles fine runs perfectly with delays executing and final calculation output printf.
The resulting calculations still wrong using float type. ????
What's going on is there a problem the compiler when using double and printf.
| Code: |
#include "..\devices\24HJ256GP210.h" // Register definitions
#include "DEFS_24HJ.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#case
#define hi(x) (*(&x+1))
#zero_ram
#zero_local_ram
#fuses NOWDT //No Watch Dog Timer
#fuses NOWRTB //Boot block not write protected
#fuses NOBSS //No boot segment
#fuses NOPROTECT //Code Protected form reads
#fuses NOWRT //Program memory not write protected
#fuses PR_PLL //Primary Oscillator with PLL
#fuses NOCKSFSM //Clock Switching disabled
#fuses NOOSCIO //OSC2 is clock output
#fuses HS //Primary oscillator HS
#fuses WINDIS //WD in non-windows mode
#fuses WPRES128
#fuses WPOSTS16
#fuses NOPUT //No Power on Reset Timer value
#fuses NOIESO //Internal external switch Over mode disabled
#fuses NORSS //No secure segment RAM
#fuses NOSSS //No secure segment
#fuses NOWRTSS //Secure segment not write protected
#fuses NORBS //No Boot Ram defined
#fuses NODEBUG //No debug mode
#fuses NOJTAG //JTAG DISABLED
//#fuses ICSP2 //ICD2 uses pin2
#use delay(clock=32000000)
#use rs232(baud=57600, xmit=PIN_F3, rcv=PIN_F2, ERRORS)
#define PI 3.14159265358979323846
#define deg2rad(x) ((x) * PI / 180.0)
#define rad2deg(x) ((x) * 180.0 / PI)
double compute_jd(unsigned 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;
// 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;
unsigned int32 y = year + 4800 - a;
int m = month + 12*a - 3;
printf("a = %d y = %d m = %d\r\n", a,y,m);
double jd = day + (153*m + 2)/5 + 365*y + y/4 - y/100 + y/400 - 32045;
jd += (utc_hour + utc_minute/60.0) / 24.0 - 0.5;
printf("jd = %f\r\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 = (double)fmod(357.52911 + 35999.05029*T, 360.0);
double M_rad = deg2rad(M);
double C = (double)deg2rad(1.914602*sin(M_rad) + 0.019993*sin(2*M_rad) + 0.000289*sin(3*M_rad));
double L0 = (double)deg2rad(fmod(280.46646 + 36000.76983*T, 360.0));
double L = L0 + C;
printf( "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
double epsilon = deg2rad(23.43929111 - 0.0130042*T);
double RA = atan2(cos(epsilon)*sin(L), cos(L));
double delta = asin(sin(epsilon)*sin(L));
double GMST = fmod(280.46061837 + 360.98564736629*(jd - 2451545.0), 360.0);
double H = (double)deg2rad(fmod(GMST + lon - (double)rad2deg(RA) + 180.0, 360.0) - 180.0);
printf( "RA->%2.5f\r\n\r",rad2deg(RA));
printf( "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
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
}
main()
{
// setup_oscillator(OSC_CRYSTAL, 40000000, 8000000);
setup_oscillator( OSC_INTERNAL, 32000000);
TRISA = 0b1111111100000000; //
TRISF = 0b1111111100001111; //
printf( "\r\nHello World\r\n");
int month, day, hour, minute;
unsigned int year;
double tz, lat, lon;
// printf("Enter year: "); scanf("%d", &year);
printf( "\r\n\r\nEnter month: "); scanf( "%d", &month); printf( "%d\r\n", month);
printf( "Enter day: "); scanf( "%d", &day); printf( "%d\r\n", day);
printf( "Enter hour: "); scanf( "%d", &hour); printf( "%d\r\n", hour);
printf( "Enter minute: "); scanf( "%d", &minute); printf( "%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;
// printf( "Calculating Sun Position for time %02d/%02d/%ld %02d:%02d:00\r\n", month, day, year, hour, minute);
double jd = compute_jd(year, month, day, hour, minute, tz);
delay_ms(5000); // THIS DELAY WORKS FOR DOUBLES AND FLOATS
double azimuth, altitude;
sun_position(jd, lat, lon, &azimuth, &altitude);
delay_ms(5000); // THIS DELAY NOT FOR DOUBLES, BUT DOES FOR FLOATS
double tilt = compute_tilt(altitude, azimuth);
delay_ms(5000); // THIS DELAY NOT FOR DOUBLES, BUT DOES FOR FLOATS
// IT APPEARS THAT NO CODE EXECUTES AFTER THIS POINT USING DOUBLE
printf( "Calculating Sun Position for time %02d/%02d/%u %02d:%02d:00\r\n", month, day, year, hour, minute);
printf( "GEO (%f, %f)\r\n", lat, lon);
printf( "\r\nSun Position:\r\n");
printf( " Azimuth: %.2f Altitude: %.2f\r\n", azimuth, altitude);
printf( " East/West Optimal Tilt Angle from Horizontal: %.2f\r\n", tilt);
delay_ms(5000); // THIS DELAY NOT FOR DOUBLES, BUT DOES FOR FLOATS
reset_cpu();
return 0;
}
|
**** NOW SHOWING OUTPUT WHEN USING double FOR DECLARATION ****
Hello World
Enter month: 6
Enter day: 21
Enter hour: 10
Enter minute: 45
a = 0 y = 6825 m = 3
jd = 2460848.15
**** THERE IS A delay_ms(5000) statement here in the code, and this code executes ****
T = 0.25470 M = 166.72253 M_rad = 2.90985 C = 0.00752 L0 = 1.57250 L = 1.58002
RA->90.57637
epsilon = 0.40903 RA = 1.58085 delta = 0.40901 GMST = 146.34 H = -0.41
**** NO CODE EXECUTED AFTER THIS POINT, NO DELAY STATEMENTS EXECUTE IT JUST RESTART FROM THE TOP. ABOVE CALCULATION APPEAR TO BE CORRECT. ****
Hello World
Enter month:
**** NOW SHOWING OUTPUT WHEN USING float FOR DECLARATION ****
Hello World
Enter month: 6
Enter day: 21
Enter hour: 10
Enter minute: 45
a = 0 y = 6825 m = 3
jd = 2460848.00
T = 0.25470 M = 166.56738 M_rad = 2.90714 C = 0.00760 L0 = 1.56979 L = 1.57739
RA->90.41230
epsilon = 0.40903 RA = 1.57799 delta = 0.40902 GMST = 89.75 H = -1.39
Calculating Sun Position for time 06/21/2025 10:45:00
GEO (43.65, -79.38)
Sun Position:
Azimuth: 78.84 Altitude: 22.91
East/West Optimal Tilt Angle from Horizontal: -66.69
Hello World
Enter month:
All code executes to the end with the printf statements printing and delay_ms(5000) statements working when float is used as a type.
What's wrong.    
Thanks again.
ARG 
I WAS LOOKING AT THE .lst THIS FILE AND CAME ACROSS THE HEADER FILE CONTAINING.
Standard template: float strtof(char * s,char *endptr)
.................... float48 strtof48(char *s,char *endptr);
.................... float64 strtod(char *s,char *endptr);
I saw above float48 & float64 type. Will these work
So I recompiled the program using float48 and everything works fine with the correct calculated values. 
I recompiled with float64, DOES NOT WORK it errors like double type. 
Somebody is going to get in trouble. 
Waisted another day. I wonder if Pic18F has similar types.  |
|
|
Ttelmah
Joined: 11 Mar 2010
Posts: 19761
|
|
| Posted: Sat Mar 29, 2025 11:29 pm |
|
|
First get rid of your definition of 'hi'. This would cause an error trap on
a 16bit PIC, since it could result in your trying to do a 16bit access to an
address offset by 8bits. Not legal......
You may be hitting a maths problem that existed about your compiler
version. There was a bug where 64bit floats caused a trap error on
certain operations. Only existed for about two compiler versions, and
it was then pointed out and fixed (I ws the one who found it).
I'me fairly sure this was around your version.
One other critical thing. Expand the stack......
The 16bit PIC's use stack for variables. Result a printf of a double or
arithmetic on these, needs hundreds of bytes of stack space. Expand it to
at least 1KB.
OK. Looking at this a bit further, a lot of comments. First, forget most of
your fuses, but in particular the ones setting the clock. Similarly get
rid of your clocl setup line. Just specify the delay as shown in the attached.
Then though modern C allows it, I have removed the continuous definitions
of the variables 'inline'. Historically C did not allow this, but critically,
doing this does not allow a debugger to stop at the line. Result makes
working out what goes wrong much harder. So I've switched to declaring
the variables first in the first routine. Then I've also changed to using
reference variables. This results in much less code complexity, when
doing a job where the original variable in the main routine wants to
be changed. I've also expanded the stack.
Then at the end you have a problem in understanding. A conventional
C program being called from a command line, can return. A standalone
'only program', can't. Unfortunately, resetting the CPU at this point
means that everything in the serial buffer will be lost. So modified
this way:
| Code: |
#include "..\devices\24HJ256GP210.h" // Register definitions
#build (stack=512) //larger stack
#use delay(INTERNAL=32MHz) //This sets up all the clock fuses for you.
//Do this before loading any code.
#include "DEFS_24HJ.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#case
//#define hi(x) (*(&x+1))
#zero_ram
//#zero_local_ram //Don't know what this is meant to do. All the RAM
//has already been zeroed.
#fuses NOWDT //No Watch Dog Timer
#fuses NOWRTB //Boot block not write protected
#fuses NOBSS //No boot segment
#fuses NOPROTECT //Code Protected form reads
#fuses NOWRT //Program memory not write protected
#fuses WINDIS //WD in non-windows mode
#fuses WPRES128
#fuses WPOSTS16
#fuses NOPUT //No Power on Reset Timer value
#fuses NOIESO //Internal external switch Over mode disabled
#fuses NORSS //No secure segment RAM
#fuses NOSSS //No secure segment
#fuses NOWRTSS //Secure segment not write protected
#fuses NORBS //No Boot Ram defined
#fuses NODEBUG //No debug mode
#fuses NOJTAG //JTAG DISABLED
//#fuses ICSP2 //ICD2 uses pin2
#use rs232(baud=57600, YART1 ERRORS)
//#define PI 3.14159265358979323846 //Already defined in math.h
#define deg2rad(x) ((x) * PI / 180.0)
#define rad2deg(x) ((x) * 180.0 / PI)
double compute_jd(unsigned 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;
// 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;
unsigned int32 y = year + 4800 - a;
int m = month + 12*a - 3;
printf("a = %d y = %d m = %d\r\n", a,y,m);
double jd = day + (153*m + 2)/5 + 365*y + y/4 - y/100 + y/400 - 32045;
jd += (utc_hour + utc_minute/60.0) / 24.0 - 0.5;
printf("jd = %f\r\n", jd);
return jd;
}
void sun_position(double jd, double lat, double lon, double & azimuth, double & altitude) {
double T,M,M_rad,C,L,L0,epsilon,RA,delta,GMST,H,lat_rad,sin_alt,cos_azm;
T = (jd - 2451545.0) / 36525.0;
// Sun's coordinates
M = (double)fmod(357.52911 + 35999.05029*T, 360.0);
M_rad = deg2rad(M);
C = (double)deg2rad(1.914602*sin(M_rad) + 0.019993*sin(2*M_rad) + 0.000289*sin(3*M_rad));
L0 = (double)deg2rad(fmod(280.46646 + 36000.76983*T, 360.0));
L = L0 + C;
printf( "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
epsilon = deg2rad(23.43929111 - 0.0130042*T);
RA = atan2(cos(epsilon)*sin(L), cos(L));
delta = asin(sin(epsilon)*sin(L));
GMST = fmod(280.46061837 + 360.98564736629*(jd - 2451545.0), 360.0);
H = (double)deg2rad(fmod(GMST + lon - (double)rad2deg(RA) + 180.0, 360.0) - 180.0);
printf( "RA->%2.5f\r\n\r",rad2deg(RA));
printf( "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
lat_rad = deg2rad(lat);
sin_alt = sin(delta)*sin(lat_rad) + cos(delta)*cos(lat_rad)*cos(H);
altitude = rad2deg(asin(sin_alt));
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
}
void main(void)
{
set_tris_a(0b1111111100000000); //
set_tris_b(0b1111111100001111); //
printf( "\r\nHello World\r\n");
int month, day, hour, minute;
unsigned int year;
double tz, lat, lon;
while (TRUE)
{
// printf("Enter year: "); scanf("%d", &year);
printf( "\r\n\r\nEnter month: "); scanf( "%d", &month); printf( "%d\r\n", month);
printf( "Enter day: "); scanf( "%d", &day); printf( "%d\r\n", day);
printf( "Enter hour: "); scanf( "%d", &hour); printf( "%d\r\n", hour);
printf( "Enter minute: "); scanf( "%d", &minute); printf( "%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;
// printf( "Calculating Sun Position for time %02d/%02d/%ld %02d:%02d:00\r\n", month, day, year, hour, minute);
double jd = compute_jd(year, month, day, hour, minute, tz);
delay_ms(5000); // THIS DELAY WORKS FOR DOUBLES AND FLOATS
double azimuth, altitude;
sun_position(jd, lat, lon, azimuth, altitude); //pass az and alt by reference
delay_ms(5000); // THIS DELAY NOT FOR DOUBLES, BUT DOES FOR FLOATS
double tilt = compute_tilt(altitude, azimuth);
delay_ms(5000); // THIS DELAY NOT FOR DOUBLES, BUT DOES FOR FLOATS
// IT APPEARS THAT NO CODE EXECUTES AFTER THIS POINT USING DOUBLE
printf( "Calculating Sun Position for time %02d/%02d/%u %02d:%02d:00\r\n", month, day, year, hour, minute);
printf( "GEO (%f, %f)\r\n", lat, lon);
printf( "\r\nSun Position:\r\n");
printf( " Azimuth: %.2f Altitude: %.2f\r\n", azimuth, altitude);
printf( " East/West Optimal Tilt Angle from Horizontal: %.2f\r\n", tilt);
delay_ms(5000); // THIS DELAY NOT FOR DOUBLES, BUT DOES FOR FLOATS
} //avoid crashing here
}
|
With those changes, the code runs and loops asking for a new set of values. |
|
|
|
|
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
