diff --git a/G4:G431/CORDIC/astro.c b/G4:G431/CORDIC/astro.c index ae4796e..c85ae45 100644 --- a/G4:G431/CORDIC/astro.c +++ b/G4:G431/CORDIC/astro.c @@ -22,11 +22,11 @@ #include "astro.h" #include "cordic.h" -static int sincosflag = 0; // math.h +static int sincosflag = 0; // 0: math.h, 1: cordic // longitude/latitude + in rad/hrs //static float longitude = 41.44143375f, latitude = 43.6535278f; -static float lat_rad = 43.6535278f * M_PIf / 180.f; -static float long_hrs = 41.44143375f / 15.f; +static float lat_rad = DEG2RAD(43.6535278f); +static float long_hrs = DEG2HOURS(41.44143375f); static void sincosf_m(float angle, float *s, float *c){ if(s) *s = sin(angle); @@ -77,13 +77,13 @@ float LST_from_unix(uint32_t t){ /* 3. Convert Hour Angle (HA) to Right Ascension (RA) and vice versa. All angles in degrees. LST is Local Sidereal Time in degrees. */ -float ha_to_ra(float ha, float lst){ - float ra = lst - ha; +float ha_to_ra(float ha, float lst_deg){ + float ra = lst_deg - ha; return normalize_degrees(ra); } -float ra_to_ha(float ra, float lst){ - float ha = lst - ra; +float ra_to_ha(float ra, float lst_deg){ + float ha = lst_deg - ra; // Hour angle is usually in range [-180,180) ha = normalize_degrees(ha); if (ha > 180.0f) ha -= 360.0f; @@ -93,8 +93,8 @@ float ra_to_ha(float ra, float lst){ /* 4. Convert Altitude-Azimuth coordinates to Equatorial (Hour Angle, Declination) and back. All angles in degrees. Azimuth is measured from North through East. */ void altaz_to_hadec(float alt_deg, float az_deg, float *ha_deg, float *dec_deg){ - float alt = alt_deg * M_PIf / 180.0f; - float az = az_deg * M_PIf / 180.0f; + float alt = DEG2RAD(alt_deg); + float az = DEG2RAD(az_deg); float sin_alt, cos_alt, sin_az, cos_az, sin_lat, cos_lat; sincosf(alt, &sin_alt, &cos_alt); @@ -110,13 +110,13 @@ void altaz_to_hadec(float alt_deg, float az_deg, float *ha_deg, float *dec_deg){ float y = -cos_alt * sin_az; float ha = atan2f(y, x); // radians - *ha_deg = ha * 180.0f / M_PIf; - *dec_deg = dec * 180.0f / M_PIf; + *ha_deg = RAD2DEG(ha); + *dec_deg = RAD2DEG(dec); } void hadec_to_altaz(float ha_deg, float dec_deg, float *alt_deg, float *az_deg){ - float ha = ha_deg * M_PIf / 180.0f; - float dec = dec_deg * M_PIf / 180.0f; + float ha = DEG2RAD(ha_deg); + float dec = DEG2RAD(dec_deg); float sin_dec, cos_dec, sin_ha, cos_ha, sin_lat, cos_lat; sincosf(dec, &sin_dec, &cos_dec); @@ -130,10 +130,10 @@ void hadec_to_altaz(float ha_deg, float dec_deg, float *alt_deg, float *az_deg){ /* Azimuth (from North through East) */ float x = sin_dec * cos_lat - cos_dec * sin_lat * cos_ha; float y = -cos_dec * sin_ha; - float az = atan2f(y, x); // radians, [-π, π] + float az = atan2f(y, x); // radians - *alt_deg = alt * 180.0f / M_PIf; - *az_deg = az * 180.0f / M_PIf; + *alt_deg = DEG2RAD(alt); + *az_deg = DEG2RAD(az); if (*az_deg < 0.0f) *az_deg += 360.0f; } @@ -150,7 +150,7 @@ void hadec_to_altaz(float ha_deg, float dec_deg, float *alt_deg, float *az_deg){ * @param refa Output: tan(Z) coefficient (radians) * @param refb Output: tan^3(Z) coefficient (radians) */ -static void refco_f32(float phpa, float tc, float rh, float wl, float *refa, float *refb) { +void refco_f32(float phpa, float tc, float rh, float wl, float *refa, float *refb) { // Restrict input parameters to safe values (clamp) float t = tc; if(t < -150.0f) t = -150.0f; @@ -195,10 +195,68 @@ static void refco_f32(float phpa, float tc, float rh, float wl, float *refa, flo if(refb) *refb = -gamma * (beta - gamma / 2.0f); } -//alt_corrected = alt_apparent + refraction -float refraction(float phpa, float tc, float rh, float Z_rad){ +#if 0 +void refco_f32(float phpa, float tc, float rh, float wl, float *refa, float *refb){ + // Restrict input parameters to safe values (clamp) + float t = tc; + if(t < -150.0f) t = -150.0f; + if(t > 200.0f) t = 200.0f; + + float p = phpa; + if(p < 0.0f) p = 0.0f; + if(p > 10000.0f) p = 10000.0f; + + float r = rh; + if(r < 0.0f) r = 0.0f; + if(r > 1.0f) r = 1.0f; + + float w = wl; + if(w < 0.1f) w = 0.1f; + if(w > 10.0f) w = 10.0f; + + // Water vapour pressure at the observer + float pw = 0.0f; + if(p > 0.0f){ + // Saturation vapour pressure (empirical formula) + float ps = powf(10.0f, (0.7859f + 0.03477f * t) / (1.0f + 0.00412f * t)) + * (1.0f + p * (4.5e-6f + 6e-10f * t * t)); + float denom = 1.0f - (1.0f - r) * ps / p; + if (denom < 1e-12f) denom = 1e-12f; + pw = r * ps / denom; + } + + // Temperature in Kelvin + float tk = t + 273.15f; + + // Refractive index minus 1 at the observer (gamma = (n - 1) at the observer) + float gamma; + // Optical/IR: wavelength-dependent formula + float wlsq = w * w; + float coef = 77.53484e-6f + (4.39108e-7f + 3.666e-9f / wlsq) / wlsq; + float num = fmaf(coef, p, -11.2684e-6f * pw); // fmaf(a,b,c) = a*b+c + gamma = num / tk; + + // Beta coefficient (from Stone, with empirical adjustments) + float beta = 4.4474e-6f * tk; + + // Refraction constants (from Green) + if(refa) *refa = gamma * (1.0f - beta); + if(refb) *refb = -gamma * (beta - gamma / 2.0f); +} +#endif + +/** + * @brief refraction - calculates refraction (z = z0 - refraction) + * @param phpa - pressure, Hpa + * @param tc - temperature, degC + * @param rh - relative humidity, 0..1 + * @param zd - zenith distance, degrees + * @return refraction, degrees + */ +float refraction(float phpa, float tc, float rh, float zd){ float A, B; refco_f32(phpa, tc, rh, 0.55, &A, &B); - float tanZ = tanf(Z_rad); - return A * tanZ + B * tanZ * tanZ * tanZ; + float tanZ = tanf(DEG2RAD(zd)); + float refr = A * tanZ + B * tanZ * tanZ * tanZ; + return RAD2DEG(refr); } diff --git a/G4:G431/CORDIC/astro.h b/G4:G431/CORDIC/astro.h index bae29a8..9490a61 100644 --- a/G4:G431/CORDIC/astro.h +++ b/G4:G431/CORDIC/astro.h @@ -20,13 +20,25 @@ #include +#ifndef M_PIf +#define M_PIf 3.141592653589793f +#endif + +#define DEG2HOURS(d) ((d) / 15.f) +#define HOURS2DEG(h) ((h) * 15.f) +#define DEG2RAD(d) ((d) * M_PIf / 180.f) +#define RAD2DEG(r) ((r) * 180.f / M_PIf) +#define DEG2ARCSEC(d) ((d) * 3600.f) +#define RAD2ARCSEC(r) ((r) * 180.f * 3600.f / M_PI) + void set_sincos(int iscordic); int get_sincos(); float MJD_from_unix(uint32_t t); //float LST_from_mjd(float mjd); float LST_from_unix(uint32_t t); -float ha_to_ra(float ha, float lst); -float ra_to_ha(float ra, float lst); +float ha_to_ra(float ha, float lst_deg); +float ra_to_ha(float ra, float lst_deg); void altaz_to_hadec(float alt_deg, float az_deg, float *ha_deg, float *dec_deg); void hadec_to_altaz(float ha_deg, float dec_deg, float *alt_deg, float *az_deg); -float refraction(float phpa, float tc, float rh, float Z_rad); +void refco_f32(float phpa, float tc, float rh, float wl, float *refa, float *refb); +float refraction(float phpa, float tc, float rh, float zd); diff --git a/G4:G431/CORDIC/commproto.cpp b/G4:G431/CORDIC/commproto.cpp index 6d86e7a..1881762 100644 --- a/G4:G431/CORDIC/commproto.cpp +++ b/G4:G431/CORDIC/commproto.cpp @@ -35,6 +35,42 @@ static int (*SEND)(const char *str) = nullptr; extern volatile uint32_t Tms; +// input coordinates (degrees) and meteo +static float phpa = 800.f, tc = 0.f, rh = 0.5f, az = 0.f, zd = 20.f, ra = 0.f, ha = 0.f, dec = 0.f; +// starting time parameters: UNIX-time value and Tms when setter called +static uint32_t unixt0 = 0, Tms0 = 0; + +// get current UNIX time by Tms +static uint32_t curUNIXt(){ + return unixt0 + (Tms - Tms0 + 500) / 1000; +} + +constexpr uint32_t hash(const char* str, uint32_t h = 0){ + return *str ? hash(str + 1, h + ((h << 7) ^ *str)) : h; +} + +// structure for setters/getters of float parameters +typedef struct { + const char* name; // command name + uint32_t hash; // its hash + float* var; // pointer to variable + float min_val; // minimal and maximal values + float max_val; +} VarEntry; + +// table for variables +static const VarEntry var_table[] = { + {"az", hash("az"), &az, -180.0f, 180.0f}, + {"dec", hash("dec"), &dec, -90.0f, 90.0f}, + {"ha", hash("ha"), &ha, -180.0f, 180.0f}, + {"ra", hash("ra"), &ra, 0.0f, 360.0f}, + {"zd", hash("zd"), &zd, 0.0f, 90.0f}, + {"humid", hash("humid"), &rh, 0.0f, 1.0f}, + {"press", hash("press"), &phpa, 0.0f, 1200.0f}, + {"temp", hash("temp"), &tc, -273.15f, 100.0f}, +}; +static const size_t var_table_size = sizeof(var_table) / sizeof(var_table[0]); + //static uint8_t curbuf[MAXSTRLEN]; // COMMAND(USART, "Read USART data or send (USART=hex)") @@ -42,15 +78,33 @@ extern volatile uint32_t Tms; #define STR_HELPER(x) #x #define STR(x) STR_HELPER(x) +// text delimeter for help function +#define DELIMETER(text) +// setters/getters for common `find_var_by_hash` and `handle_var` +#define FLOATVAR(cmd, desc) // list of all commands and handlers #define COMMAND_TABLE \ COMMAND(help, "show this help") \ - COMMAND(sets, "set sin-cos to cordic (1) or math (0)") \ - COMMAND(sincos, "calculate sin/cos for given angle in degrees") \ + DELIMETER("Test functions") \ COMMAND(testc, "test CORDIC function: sincos, sin, cos, atan, sqrt, log") \ COMMAND(testm, "test math function: sin, cos, atan, sqrt, log") \ - COMMAND(time, "show MJD and LST for given UNIX-time") \ - + COMMAND(sincos, "calculate sin/cos for given angle in degrees") \ + DELIMETER("Coordinates") \ + FLOATVAR(az, "azimuth, deg (-180..180)") \ + FLOATVAR(dec, "DEC, deg (-90..90)") \ + FLOATVAR(ha, "HA, deg (-180..180)") \ + FLOATVAR(ra, "RA, deg (0..360)") \ + FLOATVAR(zd, "zenith distance, deg (0..90)") \ + DELIMETER("Astro functions") \ + COMMAND(sets, "sin-cos is cordic (1) or math (0)") \ + COMMAND(time, "show MJD and LST for given UNIX-time (and set it as system)") \ + COMMAND(azhd, "convert alt-az to ha-dec (0, default) or vise versa (1) and change input values") \ + COMMAND(hara, "convert ha to ra (0, default) or vice versa (1) and change input values") \ + COMMAND(refr, "calculate refraction for current zd") \ + DELIMETER("Meteo conditions") \ + FLOATVAR(humid, "rel. humidity, 0..1") \ + FLOATVAR(press, "atm. pressure (hpa)") \ + FLOATVAR(temp, "temperature, degC") \ typedef struct { const char *name; @@ -63,9 +117,17 @@ COMMAND_TABLE #undef COMMAND static const CmdInfo cmdInfo[] = { // command name, description - for `help` -#define COMMAND(name, desc) { #name, desc }, +#undef DELIMETER +#undef FLOATVAR +#define DELIMETER(text) { nullptr, text }, +#define FLOATVAR(name, desc) { #name, desc }, +#define COMMAND(name, desc) { #name, desc }, COMMAND_TABLE #undef COMMAND +#undef FLOATVAR +#undef DELIMETER +#define FLOATVAR(a, b) +#define DELIMETER(t) }; static const char* errtxt[ERR_AMOUNT] = { @@ -95,7 +157,7 @@ const char *EQ = " = "; // equal sign for getters * @return setter (part after `=` without leading spaces) or NULL if none */ static char *splitargs(char *args, int32_t *parno){ - if(!args) return NULL; + if(!args) return nullptr; uint32_t U32; char *next = getnum(args, &U32); int p = -1; @@ -104,7 +166,7 @@ static char *splitargs(char *args, int32_t *parno){ next = strchr(next, '='); if(next){ if(*(++next)) next = omit_spaces(next); - if(*next == 0) next = NULL; + if(*next == 0) next = nullptr; } return next; } @@ -131,15 +193,17 @@ static bool argsvals(char *args, int32_t *parno, int32_t *parval){ static errcodes_t cmd_help(const char*, char*){ SEND(REPOURL); for(size_t i = 0; i < sizeof(cmdInfo)/sizeof(cmdInfo[0]); i++){ - SEND(cmdInfo[i].name); - SEND(" - "); + if(cmdInfo[i].name){ + SEND(cmdInfo[i].name); + SEND(" - "); + }else SEND(" "); SEND(cmdInfo[i].desc); SEND("\n"); } return ERR_AMOUNT; } static const char* parse_func_name(char *args){ - char *setter = splitargs(args, NULL); + char *setter = splitargs(args, nullptr); if(!setter) return nullptr; // remove trailing spaces char *p = setter; @@ -148,11 +212,24 @@ static const char* parse_func_name(char *args){ return setter; } +#if 0 +static errcodes_t cmd_cordic(const char *cmd, char *args){ + char *setter = splitargs(args, nullptr); + if(setter){ + uint32_t u; + if(!getint(setter, &u)) return ERR_BADVAL; + set_sincos(u); + } + int i = get_sincos(); + CMDEQ(); + +} +#endif // calculate sin/cos static errcodes_t cmd_sincos(const char *, char *args){ - char *setter = splitargs(args, NULL); + char *setter = splitargs(args, nullptr); float f; - if(!setter || !getfloat(setter, &f)) return ERR_BADVAL; + if(!setter || setter == getfloat(setter, &f)) return ERR_BADVAL; f = f/180.f * M_PIf; SEND("mathsin="); SEND(float2str(sinf(f), 7)); SEND("\nmathcos="); SEND(float2str(cosf(f), 7)); @@ -220,29 +297,105 @@ static errcodes_t cmd_testc(const char*, char *args){ return ERR_AMOUNT; } -static errcodes_t cmd_time(const char *, char *args){ - char *setter = splitargs(args, NULL); - if(!setter) return ERR_BADVAL; - uint32_t unix_time; - if(setter == getnum(setter, &unix_time)) return ERR_BADVAL; - float mjd = MJD_from_unix(unix_time); - SEND("MJD="); SEND(float2str(mjd, 7)); - SEND("\nLST="); SEND(float2str(LST_from_unix(unix_time), 7)); +static errcodes_t cmd_time(const char *cmd, char *args){ + char *setter = splitargs(args, nullptr); + if(setter){ + if(setter == getnum(setter, &unixt0)) return ERR_BADVAL; + Tms0 = Tms; + } + CMDEQ(); + uint32_t tnow = curUNIXt(); + SEND(u2str(tnow)); + float mjd = MJD_from_unix(tnow); + SEND("\nMJD="); SEND(float2str(mjd, 7)); + SEND("\nLST="); SEND(float2str(LST_from_unix(tnow), 7)); SEND("\n"); return ERR_AMOUNT; } static errcodes_t cmd_sets(const char *cmd, char *args){ int32_t val; - if(argsvals(args, NULL, &val)) set_sincos(val); + if(argsvals(args, nullptr, &val)) set_sincos(val); CMDEQ(); if(get_sincos()) SEND("CORDIC\n"); else SEND("MATH\n"); return ERR_AMOUNT; } -constexpr uint32_t hash(const char* str, uint32_t h = 0){ - return *str ? hash(str + 1, h + ((h << 7) ^ *str)) : h; +// get u32, return value or 0 in case of getter +static uint32_t getflag(char *args){ + uint32_t val = 0; + char *setter = splitargs(args, nullptr); + if(setter) getnum(setter, &val); + return val; +} + +static errcodes_t cmd_hara(const char *, char *args){ + float lst_deg = HOURS2DEG(LST_from_unix(curUNIXt())); + if(getflag(args)){ // ra to ha + ha = ra_to_ha(ra, lst_deg); + SEND("ha = "); SEND(float2str(ha, 7)); + }else{ // ha to ra + ra = ha_to_ra(ha, lst_deg); + SEND("ra = "); SEND(float2str(ra, 7)); + } + SEND("\n"); + return ERR_AMOUNT; +} + +static errcodes_t cmd_azhd(const char *, char *args){ + float alt; + if(getflag(args)){ // hd to az + hadec_to_altaz(ha, dec, &alt, &az); + zd = 90.f - alt; + SEND("az = "); SEND(float2str(az, 7)); + SEND("zd = "); SEND(float2str(zd, 7)); + }else{ // az to hd + alt = 90.f - zd; + altaz_to_hadec(alt, az, &ha, &dec); + // and recalculate ra + float lst_deg = HOURS2DEG(LST_from_unix(curUNIXt())); + ra = ha_to_ra(ha, lst_deg); + SEND("ha = "); SEND(float2str(ha, 7)); + SEND("\nra = "); SEND(float2str(ra, 7)); + SEND("\ndec = "); SEND(float2str(dec, 7)); + } + SEND("\n"); + return ERR_AMOUNT; +} + +static errcodes_t cmd_refr(const char *cmd, char *){ + float A, B; + refco_f32(phpa, tc, rh, 0.55, &A, &B); + SEND("coeffs a/b: "); SEND(float2str(A, 3)); SEND(", "); SEND(float2str(B, 3)); SEND("\n"); + float refr = refraction(phpa, tc, rh, zd); + CMDEQ(); + SEND(float2str(DEG2ARCSEC(refr), 2)); + SEND("\n"); + return ERR_AMOUNT; +} + +// setter/getter of floats +static errcodes_t handle_var(const VarEntry* entry, const char* cmd, char* setter){ + if(setter){ + float val; + if(setter == getfloat(setter, &val)) return ERR_BADPAR; + if(val < entry->min_val || val > entry->max_val) return ERR_BADVAL; + *entry->var = val; + } + CMDEQ(); + SEND(float2str(*entry->var, 7)); + SEND("\n"); + return ERR_AMOUNT; +} + +// linear search by hash for variable +static const VarEntry* find_var_by_hash(uint32_t h){ + for(size_t i = 0; i < var_table_size; ++i){ + if(var_table[i].hash == h) + return &var_table[i]; + } + return nullptr; } const char *parse_cmd(int (*sendfun)(const char *), char *str){ @@ -255,12 +408,18 @@ const char *parse_cmd(int (*sendfun)(const char *), char *str){ char *restof = (char*) str; uint32_t h = hash(command); errcodes_t ecode = ERR_AMOUNT; + const VarEntry *entry = nullptr; switch(h){ #define COMMAND(name, desc) case hash(#name): ecode = cmd_ ## name(command, restof); break; COMMAND_TABLE #undef COMMAND - default: SEND("Unknown command, try 'help'\n"); break; + default: + if((entry = find_var_by_hash(h))){ + ecode = handle_var(entry, command, splitargs(restof, nullptr)); + }else{ + SEND("Unknown command, try 'help'\n"); + } } if(ecode < ERR_AMOUNT) return errtxt[ecode]; - return NULL; + return nullptr; } diff --git a/G4:G431/CORDIC/cordic.bin b/G4:G431/CORDIC/cordic.bin index 678d43f..6dfcaff 100755 Binary files a/G4:G431/CORDIC/cordic.bin and b/G4:G431/CORDIC/cordic.bin differ diff --git a/G4:G431/CORDIC/cordic.c b/G4:G431/CORDIC/cordic.c index 359de15..97f8245 100644 --- a/G4:G431/CORDIC/cordic.c +++ b/G4:G431/CORDIC/cordic.c @@ -37,7 +37,7 @@ static int32_t float_to_q31(float x){ } // Convert Q1.31 to float -static float q31_to_float(int32_t q){ +TRUE_INLINE float q31_to_float(int32_t q){ return (float)q / Q31_BASE; } diff --git a/G4:G431/CORDIC/version.inc b/G4:G431/CORDIC/version.inc index 7630c46..6acf276 100644 --- a/G4:G431/CORDIC/version.inc +++ b/G4:G431/CORDIC/version.inc @@ -1,2 +1,2 @@ -#define BUILD_NUMBER "48" -#define BUILD_DATE "2026-08-12" +#define BUILD_NUMBER "57" +#define BUILD_DATE "2026-09-05"