From 811d8bda5dd6c8842b24b614173901f6421de5d8 Mon Sep 17 00:00:00 2001 From: Edward Emelianov Date: Wed, 12 Aug 2026 22:35:39 +0300 Subject: [PATCH] add LST --- G4:G431/CORDIC/astro.c | 204 +++++++++++++++++++++++++++++++++++ G4:G431/CORDIC/astro.h | 32 ++++++ G4:G431/CORDIC/commproto.cpp | 62 +++++++++-- G4:G431/CORDIC/cordic.bin | Bin 14936 -> 25232 bytes G4:G431/CORDIC/cordic.c | 36 ++++--- G4:G431/CORDIC/cordic.files | 2 + G4:G431/CORDIC/cordic.h | 7 ++ G4:G431/CORDIC/strfunc.c | 114 ++++++++++++++++++++ G4:G431/CORDIC/strfunc.h | 5 +- G4:G431/CORDIC/test.c | 38 +++++-- G4:G431/CORDIC/test.h | 1 + G4:G431/CORDIC/version.inc | 4 +- 12 files changed, 464 insertions(+), 41 deletions(-) create mode 100644 G4:G431/CORDIC/astro.c create mode 100644 G4:G431/CORDIC/astro.h diff --git a/G4:G431/CORDIC/astro.c b/G4:G431/CORDIC/astro.c new file mode 100644 index 0000000..ae4796e --- /dev/null +++ b/G4:G431/CORDIC/astro.c @@ -0,0 +1,204 @@ +/* + * This file is part of the cordic project. + * Copyright 2026 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include + +#include "astro.h" +#include "cordic.h" + +static int sincosflag = 0; // math.h +// 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 void sincosf_m(float angle, float *s, float *c){ + if(s) *s = sin(angle); + if(c) *c = cos(angle); +} + +static void (*sincosf)(float, float*, float*) = sincosf_m; + +void set_sincos(int iscordic){ + if(iscordic) sincosf = cordic_sincos; + else sincosf = sincosf_m; + sincosflag = iscordic; +} + +int get_sincos(){ return sincosflag; } + +/* Helper functions for angle normalization (single precision) */ +static float normalize_degrees(float angle){ + angle = fmodf(angle, 360.0f); + if (angle < 0.0f) angle += 360.0f; + return angle; +} + +static float normalize_hours(float hours){ + hours = fmodf(hours, 24.0f); + if (hours < 0.0f) hours += 24.0f; + return hours; +} + +/* 1. Compute Modified Julian Date from UNIX time (seconds since 1970-01-01 00:00:00 UTC) */ +float MJD_from_unix(uint32_t t){ + return 40587.f + (float)t / 86400.0f; +} + +float LST_from_unix(uint32_t t){ + uint32_t days = t / 86400; + uint32_t sec = t % 86400; + + float mjd_int = 40587.0f + (float)days; + + float T = (mjd_int - 51544.5f) / 36525.0f, T2 = T*T, T3 = T2*T; + float gmst0_sec = 24110.54841f + 8640184.812866f*T + 0.093104f*T2 - 6.2e-6f*T3; + float ut1_sec = (float)sec * 1.00273790935f; + float gmst_sec = gmst0_sec + ut1_sec; + float lst_hours = gmst_sec / 3600.0f + long_hrs; + return normalize_hours(lst_hours); +} + +/* 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; + return normalize_degrees(ra); +} + +float ra_to_ha(float ra, float lst){ + float ha = lst - ra; + // Hour angle is usually in range [-180,180) + ha = normalize_degrees(ha); + if (ha > 180.0f) ha -= 360.0f; + return ha; +} + +/* 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 sin_alt, cos_alt, sin_az, cos_az, sin_lat, cos_lat; + sincosf(alt, &sin_alt, &cos_alt); + sincosf(az, &sin_az, &cos_az); + sincosf(lat_rad, &sin_lat, &cos_lat); + + /* Declination */ + float sin_dec = sin_alt * sin_lat + cos_alt * cos_lat * cos_az; + float dec = asinf(sin_dec); + + /* Hour angle (using atan2f for sign determination) */ + float x = sin_alt * cos_lat - cos_alt * sin_lat * cos_az; + 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; +} + +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 sin_dec, cos_dec, sin_ha, cos_ha, sin_lat, cos_lat; + sincosf(dec, &sin_dec, &cos_dec); + sincosf(ha, &sin_ha, &cos_ha); + sincosf(lat_rad, &sin_lat, &cos_lat); + + /* Altitude */ + float sin_alt = sin_lat * sin_dec + cos_lat * cos_dec * cos_ha; + float alt = asinf(sin_alt); + + /* 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, [-π, π] + + *alt_deg = alt * 180.0f / M_PIf; + *az_deg = az * 180.0f / M_PIf; + if (*az_deg < 0.0f) *az_deg += 360.0f; +} + +/** + * @brief Calculate refraction constants A and B for the model dZ = A*tan(Z) + B*tan^3(Z) + * + * This is a single-precision adaptation of the SOFA iauRefco / ERFA eraRefco function. + * Optimized for microcontrollers without hardware double-precision support (e.g., STM32G431). + * + * @param phpa Pressure at the observer (hPa = mbar) + * @param tc Ambient temperature at the observer (degrees C) + * @param rh Relative humidity at the observer (range 0-1) + * @param wl Wavelength (micrometers). Use 0.55 for optical, >100 for radio. + * @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) { + // 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)); + pw = r * ps / (1.0f - (1.0f - r) * ps / p); + } + + // 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; + gamma = ((77.53484e-6f + (4.39108e-7f + 3.666e-9f / wlsq) / wlsq) * p + - 11.2684e-6f * pw) / 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); +} + +//alt_corrected = alt_apparent + refraction +float refraction(float phpa, float tc, float rh, float Z_rad){ + float A, B; + refco_f32(phpa, tc, rh, 0.55, &A, &B); + float tanZ = tanf(Z_rad); + return A * tanZ + B * tanZ * tanZ * tanZ; +} diff --git a/G4:G431/CORDIC/astro.h b/G4:G431/CORDIC/astro.h new file mode 100644 index 0000000..bae29a8 --- /dev/null +++ b/G4:G431/CORDIC/astro.h @@ -0,0 +1,32 @@ +/* + * This file is part of the cordic project. + * Copyright 2026 Edward V. Emelianov . + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#pragma once + +#include + +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); +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); diff --git a/G4:G431/CORDIC/commproto.cpp b/G4:G431/CORDIC/commproto.cpp index 6d071fa..6d86e7a 100644 --- a/G4:G431/CORDIC/commproto.cpp +++ b/G4:G431/CORDIC/commproto.cpp @@ -19,9 +19,12 @@ #include extern "C"{ +#include #include +#include "astro.h" #include "commproto.h" +#include "cordic.h" #include "hardware.h" #include "strfunc.h" #include "test.h" @@ -42,8 +45,11 @@ extern volatile uint32_t Tms; // 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") \ + COMMAND(testc, "test CORDIC function: sincos, sin, cos, atan, sqrt, log") \ COMMAND(testm, "test math function: sin, cos, atan, sqrt, log") \ - COMMAND(testc, "test CORDIC function: sin, cos, atan, sqrt, log") + COMMAND(time, "show MJD and LST for given UNIX-time") \ typedef struct { @@ -103,7 +109,6 @@ static char *splitargs(char *args, int32_t *parno){ return next; } -#if 0 /** * @brief argsvals - split `args` into `parno` and setter's value * @param args - rest of string after command @@ -122,7 +127,6 @@ static bool argsvals(char *args, int32_t *parno, int32_t *parval){ } return false; } -#endif static errcodes_t cmd_help(const char*, char*){ SEND(REPOURL); @@ -134,8 +138,8 @@ static errcodes_t cmd_help(const char*, char*){ return ERR_AMOUNT; } -static const char* parse_func_name(char *args, int32_t *parno){ - char *setter = splitargs(args, parno); +static const char* parse_func_name(char *args){ + char *setter = splitargs(args, NULL); if(!setter) return nullptr; // remove trailing spaces char *p = setter; @@ -144,10 +148,25 @@ static const char* parse_func_name(char *args, int32_t *parno){ return setter; } +// calculate sin/cos +static errcodes_t cmd_sincos(const char *, char *args){ + char *setter = splitargs(args, NULL); + float f; + if(!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)); + float s, c; + cordic_sincos(f, &s, &c); + SEND("\ncordicsin="); SEND(float2str(s, 7)); + SEND("\ncordiccos="); SEND(float2str(c, 7)); + SEND("\n"); + return ERR_AMOUNT; +} + // test math function static errcodes_t cmd_testm(const char*, char *args){ - int32_t parno; - const char *fname = parse_func_name(args, &parno); + const char *fname = parse_func_name(args); if(!fname) return ERR_BADPAR; uint32_t elapsed = 0; bool ok = true; @@ -174,12 +193,14 @@ static errcodes_t cmd_testm(const char*, char *args){ // test CORDIC function static errcodes_t cmd_testc(const char*, char *args){ - int32_t parno; - const char *fname = parse_func_name(args, &parno); + const char *fname = parse_func_name(args); if(!fname) return ERR_BADPAR; uint32_t elapsed = 0; bool ok = true; - if(strcmp(fname, "sin") == 0){ + if(strcmp(fname, "sincos") == 0){ + elapsed = test_cordic_sincos(); + } + else if(strcmp(fname, "sin") == 0){ elapsed = test_cordic_sin(); }else if(strcmp(fname, "cos") == 0){ elapsed = test_cordic_cos(); @@ -199,6 +220,27 @@ 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)); + 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); + 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; } diff --git a/G4:G431/CORDIC/cordic.bin b/G4:G431/CORDIC/cordic.bin index eea0d9593609161a04084af87ced79b6215347d9..678d43f4f48709e9f00fc82272c23d952f3d97d4 100755 GIT binary patch literal 25232 zcmdVC30PCt);GS-Ndh56KoLbm3`7ON0mTsOj@DHr`u;L$%u8 z>(Ifewbk0m)}iWct);QrRy*CBIJ7EKxV4&^nDj)j1W@w*)=8k)_P+1?z0ddjpXbl> zEcQ8juf6u#YpuQZ-fQm@W@ehd1~Dc3U;I__OxaiVe*to)bp1&kYWn|*&VSS9w|~;c zpV4paGL^OGq5boJSFg3amDYci{wLhl_Wvh!TWPef^ZRFduW=1B4Jt6Z6HO6^W0p0z zndz|GZt}Pzn)+6Mb+BJSzx)_Ozk=laV#8)*(!pT&hV&+Pt?9UnnR*v&H}-AI`4l-1 zdgTRW;}QyY&DV=ce%Hg6x}v{#&yOuE5;^rsFFW4ZaMA8-*TVf4(LK(jKN@vc4K5+3 z>;qtRDCmH;>)od$Za%(bL0;?IYRdulT+=7uB)FifJ0g9by(&CU_by`^rsj3`KuzWb zUpb55b!I&B`Nc*qt=LdDj#&stQBXB+yNo|CFQB&gf+Ep3IA$hFM{qjGRufvV;y`yJ zx0&0%Er>rK0<5WJ-1hE4^7Fy@k)|w>lSMA(3tml_YTyqLO~!bs9?@>GGmG^)KWM#v zfHmSsbhX~NS>NfB+{k6foxEC`t@ijvakaB@l-5^-NG~@~mF%^B=HXabc0?uX+p(Hy z8F0C5SvEJ4t6i^GT4USrjss5vJn)<9|6ew)>_~(=qEsyy4l0CHK_MnKX5NSh<;9-O z_aMc4Q<7o&K`My}l>D>FUlP-#F_cw~I1*kMp09n1#oG*z9F>Va1!Nfi*>yR`X@4}n zYScd7SCESnTY5hpTo9c1R(jHo1eNpgn_KE^XUqAFM!kY4v- z#krsYCFCSK0NZHxHJtXzBnk46pvO*fqq_aCd$S<3p)c#Z3dbyZ)K$BUGmR

-Eb=Jpog{cG zNgn5}hiEtJk6GeNpNZKXGcB5>h}P%5{ZdgD#8Mo#u%fVBP7Ap#9x1FszOUgWA4|fM z;c|6IK}i0#_=_dF_?#?;kWUZc z@t7lejqxo==zyTi)(aycp(pMf(Fav)jC-6)KG><`LMyObm) zR$PJ9t{8r`lMLuxA&7R^+|XZO+u(7W0`wldBNudsx)gPp@V!^Qp-5!P{9oTmR@U~4 zF=s2Xqpi!0yPY92c}nDk@Z2xF?9~LJui$I=n=dQYV`WKBo3A0?M>SjVoELS})*3VP zNvUL0sZ2-m>69N?uvxD>!=;~zNf2dxM4o(bGt{~mFf;wN#p+bTjSXJxN*uOi3K&TMBpTrBAp6I$8$c?e56x%H#2Q_HTljD z;ah85y-dCGhr9cj_~uC2aL}xA#hXsKgGz(p$v%X?o8x}dl<4kl+Ubfdo#TFE>J#qW zrrz$S@w}(>3HQ?|UsYP?R;r$GUo$=7UIU(92TuwWFV1f~;y$h-@11P{hAGvZs8Uu2 z=L|SOe6sh##(nUBd2F`_KNhxxScfrIlGDwWsx#l8H_m(y$Bup{D=2#JlhG0ROuNF{ zeQ?szRHQM3`I|ACvT+fILkndBf3v&h^=4V8Zo*lm%(+s{c~bdZ&O$zb#8gd)zu^LWID2OH_-1qj>({>z>pwg>$v+ z6HiFE+?gc{xo||!UvS%o_PcR$kyxG;YA^_8l`Hfe3{BRuN`=4@LTcaluo>c+D75_s z`jB>zTKnn!;C?dai6HUfTY8n@Cu7_*Am|${dRjYgX|Yld)=@P@8z(gPcnVVE9l*n*u4jnz=5ik zB^NpL{JyF|DO6Kw(PI>Mm-yiNi5FLcgWcekdTtPQ+9e5$+q!jYH}QyU*~$Y7Axw}9 z-H&#v?Tp&XCHeJF3AJ%UxVJSJYlXg5dW_oPDqE5uJ0Mp2b2Ft}F1X&xZD%Q-s^Df> zT^jteSjjAh-8#v3%$_Ue8rNS;2i(W88kj51zpQShxz)DBBh0f07?06fL8_svd&kUr z*SN||OpCRIN|x4I(r}q2K9WCw!qd&N*vnnJ37%soo$g*8=Ia=h5WJx8~+u$HL>A`}BijpkBny&5*`7{P)PP6Q+E4p;EOJgd- zJ2&D9_Zy}eP3D%l7!8XYVL`Vs&bnj+KRrU|CcN1!J6WYn9Ypr}nzN9DrFMH-f77L5~$Wd-^gdkR8zD>`_ zXDOS0ma+7L=u0QPuk$OtZ*>}xA0d>DOZD^K|7R=ph46pHx`oCe^@V!ym0fmoPeN-g zs0ex@3pr|ynl6IWhUYDInR3UbVZq10_AGpjUt6sdGGSb z_*>^rdSo&#QE%^5&^aG1B^Z!F2qJG-K^}c*rkv5Nq%u6@=^ZTgv<%73^uspfL!Pxq#S)!vOpTq6=qF}_hSE) z`)xCq?SlL+z8}(8j#Uf0`GQ$4aD8fA_nL0H<)yILf#mH7hdka~=65qzPAvUD%7K>cSyukXr zZZC}02;AGghOCxV`hQE~CbpW(04F9cV`?nfT$7FUKX20rL?@xG=e6ibCtHoHn~C*# z(XJ7`w5x?myUfGMKCw!*)3}sVQtf}(S^pYqjVlcN{KGyVS#4$g593{Bk_&tB9)kCG z@g8J)(~j8)xm=$Gc;fNI`WT@*hCfG(UMUZq`6fPq^op{(XZpJ*WTK)~Vn$ zsIho%BbUJ9AD9_qIDCUsapEVwJ^JTdc-sD}Eu?@Y#2D@`_!C@`vtO=2da*&*Z2H@TMu%t#{kH{3%~Y0`ozKa~}AvclF=8I1{(sV^iJ5At$cB zBf@LB?WA|cPV0$R((wabKd`m(>eu((GK?+i@n7)CrTu1+eomj$yZ)mDBk%kt2{Oxn z&)58lf;|72DZrR1C@m-~C@P##IDW+e?v`o3QjiBEh!GPeb3TL-qqs|w{t;6PX&Q2V zwvF{mQ^OM225!UDnb>6-YFtCHE=(^?&&?R0kw2R(rPO*?!&m=dr>zxer~R?LtKlg8&F{xJ>HBx%JeBs~(t2AR*pp)~ zPJ43l0e*X!Tpi@1hh`ZF0{Gu@-5IVsX`?ouoEZ?O!&giqV=(MB3R`(MbIVjchglJCU-=klHS zKalUl|7Y^etI?NnSrfAUV7pc7yM2X!Mj&6YD(8l6#P4htexlUxEq3B* z==h=eu-(Gp6OwA%$G?)#r7ZiSeZikd5wghv8T76qpV4Z$&q|g%7qh$Hdc&>cds^#| zReoD7V&13HXLGpizrfeFpViy9A9S%b{zg9G8@pEcDqydza~1i6yXULpSX;8x_HjD$ z4}Y}T9V=eUC0|_Qda?Ay!WSpJ1c@?uVNSqY_7ptv%stKgl4+q~cX?p{4DoD^+xBeR zTt=%EE`jryGizLHu{&P@AG2oaP5{z;0CwqX5M$(S$p;OoZ7+$O>ekWu(E0L$m}RFR zHI`gvTPH;}*GXR4gthbiy?T_NKGCdOVSb}J2tJGIP>TkeDuL($(RFe70^KV_}=~yI7mO*zU|f+=>r86 z5xqZwD8(2)uDZA2xP5ZL{g`j4H%xofUhlptwG(%Kyws09rgg{KQ_2nW;s88WXTx@T zt6uf__xB>*b~aG|UI3lX19UPZI;($2XIOyFPS-Hdd7_O@Y8xGD;Ze}}DL|*MMCWyr z{w{ucwegb;I>l{t;@ar!bj6v*+eP;*NusM@spa+e?hm`mR$F_~v8~ntceZJ8tBwBo z8ik@w0@-gPA2pFN*z3fTQap%a%2GTTd!ygxuB9uU35@+`N}nlwX2QC*k^e_Ky?ZV5 zsxWs?Jad~J4+@1@U*m=pqP^1M=)`{ znY=`+wSIcx0kok%)}H#k0qO_bvq3%Q=8tJ-?Gu5OezZMStHtI^7JFa7Vi#8B7Ery9 z?AZnP!D{K5D7@d!PDTgpWZhk@j%e%CNVNKVo1F~Kzb|07nEg$E1D($T_Ld^iSp|D* z&sU!|I!T~Y)<#Fuj!wL(w9VegGnUsswm^dEGFA)T{lfpS4*}xMg3qVo-@h*s8U*9|Lf-5P@`5fFIjyc)DUXi!3Ccm zSiseC88R!gsEnJq#&gMep~kFwa~hRzkp-ViY*+r0tnqw8zBX_3f-KQoym>Zg)gPY~ zYABdDO%Z%Pwhc3X-V{aSxvqIJ1z&>-aF-CZSGeGFk!?6C;|8wre5d@dyhyZHJ1a|+ z%OVZyfi07*N4lu4Lu6Kc=i>DSn|DhpouGV*l%974cRI%vgx1o!rNtM{SQ*nrPGN|Q z^E#ZFpeMl6Yzp?CKah8zVBLZPg$Hx>h_$bq9b-%oYxx+%R7z(~l~Ua#y=b1~f2O<` zqh2%sXHWj8-nC4jO_=qyG$MB8Y9HB^ziWJy;nilD4AE{~mU40hqYGxH1 zw)Hu#-E3$nSzWMf=3AN%4H+fT+EZROq8t9pM`YjZ0C{$t+0u~M{8WML1r5$rit5?S z82mR+iZOh9olSXi!d?Mq6)RanR?cwK;Kj!*38lLb?ceQG%C<)z6;9RZf+I<{WGlP{rSt|1FX$E3E&W+tb_Gk0sg4p9Xg?>P-+btn zxtzzT1-mJ$(B~!#Ee-RUSL^xhOiZj`Ns4M+wl9gh`06clnP`|>ml@Jcb5;<1Gz6B{_U5 z$>E9-a!}yZL~+{Uie}oE_1edEb%)L6{ok^b_m*lys(Krz?UC0VwTZWvn%>eUK0d-W zG&?#kI(M~KvzUqFZ}*pM9HUo<;j2+hM2w^*UUhc-%m}7rxxN}-ozLgo%WxLq-!t>N zB%BF-erVT5IveUmXOsVCwd5D*OoYzba2T!mhG&`(oXaZ@wloZCUSee1w#OX0)l{xN zd&%25HAY^Q#n+p8nW&RRW_{(2P6fv1OWufyQT1#}N9!>U3pZJp;H(oTNTD^?py7#2 zjNf@<6Sl{Mq>TMlwLi>GI0Fd>l#^1*fBkmf)B@SI9)=mEONmZ~iblzrb+h_Z7{u&>aa z*YffF&!!}M9@4Lx;_ba3mn6)yl#P(^Qv#Y}V8Cn|dYcl7^ZE4X}%T*~#OYbt4RNf%&^CgY_)v?cUEe&mK zX@W;YEVkEy7{d~yZ1a=$v~0pIuRpgB=V4{mYb^~68~rvi?dVQd%yLe9?mEx+6;`?A z%e9a5;z?ic$GODMu4!`P$7cw8nu8O4N*T@*f}Lx)UCvBCB*lL&@{Bz?H+uY@W~Ejy z6gkc+mCm&sFGZL$Wg*GH87*FzXXk~#USpR1b=b$X{vl>|8FAPffG?oydZf1j#{m?d zB^w}H_y*6XP)_k$^O2_yu>H#nBMr-qEQPZ%?eN-Qk6PAESY~W#&^0eBfp1F7C^5}j zZrrAqTbqr(6z)Q^F{2~{_`Eg4m_9GYFgkwOJYF2?TUN4s9!rU`E}zGxEyr_rjOor6 zPRq1a@3Rz{b=ka@1~q8-y+&Wd0iU;|859XK z5`1+AUtZ97zWKwFkLESc+oli3=^^n!{9$#j#fZ>G3z`7ERe`=m6?8X7NTFMGOe-p zC70YpJHu0m_v~~Hk@gd(kT)288H~Q97Q`CA$2-<`DkgevQ#>wMLvJk&ottGswD9V< zMBpX?xA$GR+wa0{Y4YO|tzZvMieWV!QJW~;X}I3f1vH{Tqq9Uq*4D<4fqG#m?}&1> zRL-}Rp9z!)p`2t-@{WbPH{cw6NEYo0N$#{P*3z*2_Og=xwSA=A>utGTAUB~lUdsKUEw|&HWhGr}dq}xw z+j1u(S6v$|<(_WKeHpnNT5Vq1{Fjwj!;t-xej;ka(J!*L4(K;&@u5b)97F!l-$L|h zB49FLI^bad_1%c{kNrK-@>_omZT(G^`1;FTe2u8>C*?lcmOG(^q>w1(&S}f7LvBK? zR?3~$mb(R(*tIrJ%AFX<^)(Flk`!X3+?=*tbK`%KLIipn0V#w9q)_GcOF@PF(7Uy6 z*?=>p|G6yu<7G4CPrit=$+V+eJsPC_y_Xl3&)ck5w|OPfN=wtd?HPi3LD?}zU< zqwTc-^ASW!0)C%7Q3v?Z`(rbVRh-G#Y1vd=hA}d`8%~Cu)J|;%cu@v?>bgL^`|_6? zBMlf=GgY6ZcDrc2>3E-wvBaaEI#SY8^LO@59C zMJ_GFKp5l&DTNojm|<)#0iKL*X80Q3^L<#N)ULrQ`q;ct`J45H)?@S1kxsvGY~Dzu z>3AQGbTZzDAx*;j5Tu28PeVEe?}LzL<2@Cr9Pa~>s;n8X2&JqMGmo#~2D~V>F|1wQ zVDy^u`nAg&2CJsLgm!sjTKsv@?eabj;DxozI}pH=*U~5_`(c)sKr1*!ae-2G{v`bO z&s|AvbIW`CK3uaK@b@JtlL!9&X4u=yj2*^)YYp1~{f&Sw>26kBe!R%CHtXXqvLmtR z+v)?#i@Xpj@ZyErbh}F--KH`MOdNHai896|z3gJ4FT42eovnOQu=UlWOiXKXV)v=) z&cp#9F|l7$@r{TSnTKUDobvhWGx34uIynQFO^j(cBWZhx6QF`O8l{fhH%1U0m*0k3BmWtlu zwK9Ve_Zz5>emVz0hlyFf!vPxM-{GWPs;|0?Q!nXBY2(Guni8}-T_Gx&p_S^fCVG+y z@d+^oS#IUeN*xoI93|RjUoCXlp8kTyR{E_x{_N-RTr0Kb|CGlObf<4>Z;@pm_Pkak zd7`t2J-(tT7g+x%s&~33s76$dFW7-GbJ&wTWBgfK%_X|uM7$^u(ak)p6!H)^%FCaF zJ4I7t-*B%sU32e26iGm?F!fD$pQ%6LR?$vZHZZgE_h3zN)%^y_^(fa%a- zyb4@8W7>nrP$Hs28P<%vQTd#6C1-Q7(&2~}^UjsHi$TwGO|1Vg++sLz=^)MxM^vV| z8U9lg2KS|;v-`t=la@ECac@w3(VK!kL|xc@%PQNU?{h&!#FApnueuWvxuhEtr|`B3 zwxb1|;h;RCG9EQTFErq06U8?PlLzKp)p4}&SLu!e)sbnRus0w@2_R9$X(@LB?udQ~ z+pZb&{k+wq9=KF6I(}DFlSHDe%!`^N=4)H@W&kfAC%h0y-YC~7ih4z zmfo+o@vsxUptdPKbrAQ==o#mDMfcr^cqCC|xgUC-oH6Z2$Gxi$JbES-w+?yn_-&TN z#Cr2!Ap^HlDJQS-$=zqKGcDa5;YEGp6Ve+mE@(;+QjjyU@v%V8%s@_RA|!yTcf1(|TXLvp8!x8gY4nA+Zm~#~ zFf@;8HV9jRaSHwBxzx6PZ}jSg0^}e2UH?~m|JeWK?fSpF@wfiJEIn6)E-$_$-6B@d zJ&(Wdbah%@<9&F;b1r^4_V3D*7^^fpQQ5f>Yd^W(^j76(2M<@i;*KBl-lQGoqNg^E z*%0-0w^zBQF;Cp{X_xt*-aEE%3jOVoY@9KBAQh3;0>rU{)6OEj4mbml?zj|9xe-wS z@0hjvz#&*1&c3^^5Q5J5b+x9MEbdt|N!gzh;SCbP1=5xOUFpq{q$di>2u;aA0)Pyq)(+2 zg$F?QE8lMxcGK%`>p1Y3mKPn~Scx8S>&HGMVf)8TwWJ$=56*jkvq7Rq7@WAnrw;J* zxfd;KIdQYkZ)=uz_&JcBUqk*lL!Ko6&mr5|v41)Luc_%5ua0)^NLhF7>eMlRZ8><~ zl^&I`hJSs5fqQ5%xHP9*Dr0-2f|J+r*2{OV&!G_}lAB$aReOF-ae4-!sqz zFFp>O742|T!0G3mFX1$`e0yEZb?}T($$FDT142yTfcV=haHDy#H)4agY%>eE&AgUv z9t__RFDTz);;}nyV36?6r4nO;k;6Lc-$rn+-6cJv4_ZOL!wydyzX~fg>Jw}06Um7a zoj%G%9nNF-TN)?+J79ApGfuqZYcDBIJmVWCNs1G{_05B~&TV|)M#TCwzr@jd=)Yy8 z)GLdR7n>V{(Dxvz->`gbe4uVcORGh4;@iGBDR+giUZBy=iHoIDvc_+(%VAMld_4yc z7kQVp1&iZ~y$Pv%^;uSu4(ZBowNJyI33r8$4WvF%ip~c^%<4$B)dUH7X7Z zGFMH<{u#C_NBU)*s`#WcsyHI0$kBVsq3V03-uQdcp`}%NelF9bvE}bskam->9rEhy zqn>bLZ(qJtlWO>DyYKi*qJ3C}LmCqU@`&(#iZ)#_gZdM0d4p-cUw8gmXYaJOS$#bS z(rj!Y$1oGc+374My7Pj|y^oZy@(v6)O08*TQ4k~@q zP3QHgCcAqz)}?#!{xqJ~l;TXi)PNPpKTyw(8~Iobs=jxx!7X?LRwaIFJwUAosPzDa z9-!8vG%la$ZPxE}b;rs!aq1qdncu(~((X>g-TF66&tL3{T~6l^w#P7*JGQ+H;5Q2VMwOzcczJYH8 zQjgZSBB&H^tI=wrJI4h6N?*l_9d|}czsH&vEBAmjnC0x30|pKm`2G9ce;F~1>oD=? zGi#Njz}F{$F=dDu;H+20ZP7qaQ{!kQn{ftQ#z1e_fB1@qoY|otQN>}Ud%v}#=qqBP zPZ(n>a>)BF691|wauljR#>}Vi_=b8fvjKB#(LD0QbS%u{8`FF~%JWqC6PJ*~GP}#0 zIL{tunBrAun0&X+7`Nx>*Z6kjv8sO7L$_G?(5lQ@tdmU!Xk?B1O%>gx?P&#ljvJA1 z_BGs9rFL#z%45IUuSPw6s4D7sT^Htx`LwPJx7$-5mRZj@T|JW-<5+*YjVPy9h1wDe zWX?%THoKTa@9UUcODOKl*V(5@{877&RdL|4JMk!g9h_-1Bc)aLC-(BuTkKU)M4MVd zxkNv&JNm$=FWiVj_E2H>EhY;CPGnY=56Dg!skq&g1-Urr@D$GFYhjUMYmNJ{wGk%5m)Ff3_ zRTSDuF3*^eKP#K_gsW+{qjvB-?S*7cYjMJ{*1MZCZFI`F!MxKo1(t&Uu#jTNHxX}T zh!XE~72$0S-fHkfvR@=dtLf(aNjvS}%Er;YtPJ@yBFyf@d?EMKTlaNRX7jfa7V}vw zd8aEYy~dS=8wy!d0qJv7$w9Iv{GG!d2U1h2?(#Q+I#!a>OVo-xL5ggH%Sae4e_vbNQlF^dg z^d{~x4#!^VC5%(*;as$_%5>!@Nrw8b0j8$(N3=_{33HQ#>8BH|+Ubf8+^jUK;4LJ% zM3(M!MVfxZ?sYhP&L$~W)0~fifGNoT_YIQ5)an27y|lLLnBHmWV|wRp8dGGA9n<%3 zv#^&;9*c7sX-{cu=sMi|=%BIEcYEpU^Rar$eqR4HzUFxJ;@bS@^>IS43$&|bDK)OA zO;6_~h~Ye<6((KXHS;<>^(Vx#(^UcKR-h*p_#WsKxVi?Joidqm3zWVII!^a=H9;6Z zCQvNCNTWFK%#Tq|r@!tzQA|CpgrsU*OH1h0KwT!-?svJZ>W&&o;{5>Wg*VU+%C0Coa7JKF))OgU}V-v0sh%&+)F)AGWgu zys7jCChpy7 zQwpt`;x1?BXU6rx_?hI5*Jwqb3_jdSZIb1%0i^lW@Nkx2AKyRO zSJm~T^EFxUW-e_v;#O7cY;^Fp(9MK*E=IoX`{>wn-0~bAF=? zaefoh;a2T(D%&K?%BFn)OE%x!Qlav2s%<$=&4>)LjaocMktCyLF>uouYl%0owD|HX;#G+J{ghEZkXUQz6%y<#dx z^~Dz2akGRXvAoz%dycfDWa9n4`wD24NBiCB?P#xR?0|36IPrM%3jF`NEtbNGPc$-V z2cfHb-psbE-+3VJ3s1C+My|d|CbA}3-IvkPBh~-y_O`- ziL<0VzIFwJv~cQEt$v51kSw9tEAMd(7RT|8HkPatBP!%*fy>ZWIm(Dwn z^PF>?ifx+em4|cUj3(SbN%hHuM5JY2}uoQu`f*A;Ww=1z91{foT}-=!kiZC*kwYB9g(k5lLW8m&>wBpUiopwE90((Zm;&T1G_!ZJod!->b;>J-;mHv z#ok`2vyO&(%voH9`fWtJWn@|0sOWUox6vB933F6v=oZ}Rxx@NRyAiSX_|0!C@Rgdf z$Z?v(Xnaz=x3CB^qI0c+^$BsFR!;HCJ8;;a)6S2Hr&}@vBX(kgr6>;*kGH7zzO1L6 zp~J?uzUEP?a`%3f{q1!{CPg1&gy7A^!dF7IH)6wX;V4$X$8YX7Lk=l5UaPq$zWyS+ z*)BxXHvG(My^`gPliQivbwj;EP43l2;1~%@>e~CrLb%!ZpD{WjW)7-oj9;Uwf!U_7P9R! zF_?9L7pBPH-4@kFvZ(XdSo|Wft1Nh9&ehg^KP{Yoi}R7ar^5CT8-C`&X!I&^oA4M~ zq@64iAM(AELmuoFS|+=r)~}c|5Y=Pi1fRsk-1!Id_i?>3_O#qkCg%9w)q5}T+anKB z%%Ph=IQ;K9h2%{3&{xO?%Jg}No{$wVvE(-SyR{f|es7o8NHN7-=2d|lCSJS4iIWKwrcdnW*q_Vpe0Qhm))KMi0yKrm)Z#v%I~^e!RA8n=HA!e(rwO zkBMKKWzNqOSSu}`bfatcYYaM_j3~?Ow*<8ig;*{;^IU45qJ*H|vI`L?^34>}^Gn6n z;kdS{%fjj!O&2#p^L$<&p9pd81N7YcE9GnxzCHyjOFkLRM ze`g=r^3H4ISGw@LG`p17pLzF#tGl+#aO-Q^zTCZp@sge2UVS~6OFO9ALGzvzz0zL1 zMsG}uSc>*L+|%;=XY$kHSXzhlcgfdR9f%d47i61N!bb(ht+5_?IMyR&D+Ia6H}d$+ zdslpnFLQG(YL8r26~SX2Qsg)p_O3Hbrr6AH@9GR;^s5P7&T*+bU3*RQc0^~#z0Fcy zw@mhGCt?*fmA)N)*FyUo`sI|BdiqAwZDRwE-~4ffY-`+G%vr&T92Lr^a8pruuTW4q zD>(8Q=9{O@+>pci9IrfsdE~6#*1?k$8tPml`}TU!0C~nv*QRYXt_|xHTjZNnPNiJf zhcYF5>07Vncn4syfty`|s>rcM5mmjU+f>2e4l&WLgjVKDT%5G?I;4No{WGm!Rh3x3 z?gD2{8@IvNaJ;D=`!!&%(_qEN#1qXo@&2lN?mBJ@OSsZW@7iUEK=vlX1u4 ziKcnuYbP|xZmY2#*oSZKj(dXnZdjpCXyR{T{gfPn=K|ocy8W~tkxo^YVlLaMS5-fl z&9v~iJ#uog9Gna*0r(2O0^h}DO!I0r%VT#4v^Hok%eQgySOq-D4PyP9tRu7emu1@& zXe-QPo-praie98=E1Z7|i+eA~xiVPcT%q^@I~9^WNi@tO54(y!*5ITxq$*h`7V2(r zrGtU7TTeUM2AkT$sh`cU1Lw!;aABIA%vu>gY;(%&shkI$CQ_T*zb5(37+vIBqX{vF>mExpTMSC_n zS->8~bo}OyikSo>m&DvY$E%4 zwUOHrvTX#k9Afv=tGxV*@V$LGC}=fjI53P zd<_pa!e?(In?HUtp~4U=yf?#@QCh)5@mmB6rj8yD?# z!^{rZ~H7dy!+MGWBYhST5s}39r2V@9%8~F3nC+zzR0KZ7&u` zAKeeqycB6D!Y_kF8YUn-pkqL@^LQ;bQW!mxDg=Lh{wuf5zQTX(Z%@r(L6M0AH-?&295KB8(-U0m_7VOfr# zDH@?)fuF`q=`^^kQg_q$=sJCkEkQbCW#ZG)DLH<-q<`1fY*--%BaBL?{tV~;baq)* z87bxa&movNM#`Tgqn~0Q+>iJH(S$TXQ>AqT2VJnY%IRtVBMKPb z0wbh~^%)|Z#trdMtlYknB=AprFGLUfWOcB?60$$CckUmA-QN%em++1*GW!+%`HS3; zvvzqQ>+`O?RZc$o#)gR&iaDP^)T=zBOG*M_L)~GC5kg#5OkTH4m0pYAn&>7)eOtf$ zim7JFa`c?>*r(S;XKOIG`tPDpL@5+^Uno{w<5J&Mi(N)cCz-FpcTgzph#T}lNO8+Y zbrbfA-y2~R(;>a!&-~^4bYB#|ik8yp;5w`-QMYfj1@|k*F~$`{8d3rY_eQ%ad!wEqu0uHyM}S>`{GyGZ8c`{e+F;4e zw&)mc1pS5@G^8%9ubm(3=s3BDCDsu+<#~rX`~?6LfAIFpu5ieL;vI$eXdQ*c6CK>p zXv-vrY-nb^xnmY=VTAS7X8xbhf#ToI{M8^tZ)=-nSMeV7?`GM(J46;9ZJ8i)VXu39 z@_*c7&jE)$chqYwAK$87 z6nEhh4=YkyW;iB>O+@S~(<#5FtP-sVf`j994xWKK?yz&Z{DkG4Q?1B!P&x4i4#K!m ztZuKlvu9ZN50E4xLYZZipLv*EGivB1^^oqD<(|&+&pa9hYf{u{n&cwaq_{~mfvJxN zc7*e6S$vcHSG_Zcd&rR$rnG1rXUm9&I#5rJdJ5E&x7L&DWO~@|@TE3xPH^yHXs3?a z$#U>A&=!d+q&GEam!ZA(?K}_Zh1-B3=+JYH>sW8@skCsT8myU)AnpZ8w$M4sC$TK2 zOcv=cAKhSuACZ-fBR#62&u9y2Q*E8-$P6OwWd+%6BsXdsa$;qbR1##N(|9WFW~G|E z4cZDxV)y))JGBfMTaIQ);{BO&+vdRk3V$quHT;^5Wgj)$k0J7!O3;*aJc$K=W9oKvQJXjy=`R=;dB&UC<~ zY=widmC$zU`>_6NEvc0qMOIYg;06`d%o}wd-%ou$1_YM{i<51Kjb11uoFKs`j^cSa1%1 z-fy)<0S;$4HclT~e_>w{e(7(9t&xvn#PD zhbJEugmFi|IZES?((zvTw-|vKY4fLJ6oMi>b+>W{Fy_^2@&;;>8!wGanr|Yo!=s%0 zEnCsQ4f=kT3tbmEu9WFSj+KpTr6qZXQNcJ1&yr|=M;>C-Su1`6FO=r0RnYZJ;aPZ; zS-5Y{EZh-a!w2|n$XpB2=PJj0l`daHy$|0lXqe_fud*xNE!HXOOvm$!>X$7!*yQ#v2oB) zs=|g6ozk)h|6B7RPtbecQ=8Sxg0u73%bu3$9WFxGg+-3ze4%4v8Gco-NatM1Fj^&9 zbrj7-@wQ0Mh{YdUfCoyMG&4CSO*(<7D~;5MCN*XdseRNS&hUMlR)(>wbEp|OtAl)H z-(XzLEDk1roE8VIw&-;AI&exfb#*D=hC6j&!N&2Rnep3LyBS7DF z({ExKXU{JqB1?>mp3p2c&Re4K=Pa4MbcvKdd#Ps0yoLQ|En1>kx=1r?(c*{a&C(Z@6_ShJvHsZo=Z+L!9jTL=`!W|hpJ_1OH9rL(ElK&mwdh;!#XK6|02WZ~TTvo-S; zY95|Fck%4m68)vKmn@xy`jl$2iY8^}WohO-ws6+cd5adN6Ah5)M`=H_ji-{OC3sx& z=;EdQH1ikDmFg{MtuOI&SMC2$3+>HYFk9*$^?Q7NHkf-@Q&>E?9lKK|v;%bYkWtDZMwdP{r<>R5W2+Va^27OIGHD$&;o?k4!(M_yNi*x;JN%P*`*y zJ!?j5m~rXSvL)#Q2hN?h)cDxU0kak@7&!akho79iVBnIa3(^KJDOpf9fA*4rOCh6y z3rdzOoxOPAqTOG^C?zpKVt$*U7f6m? zRDSJAMS*9j1Q`qujARQmM$)B}F9`J02N=vZ1=FzYjA)GKn|fh1@W$#2pV9$0Q#U_C z4}87u*!9GfSj9qn}VxyOPSi{I*wrSEj*8GTr z^*lI(DIa`}rONJQ%Wq@*w)YWyWcn{QtbQB2@^%mQP{~H7#0kLat+UyyooBG;bVg>r z;A4qP#AO+3Ka)N6G0MNpWUrnF&X1Ywgd2FxnfSV1&(c+TwlYEw zI(pW*hn`K>>e-&YdeGFfVZ-$7sZn}%O0Nff{Vp)k2v`XqL!JPjA4rJ;kdL|zNCeym zumNrXh5)Jo#{pGhESKLCyaJ_nQlJ^=IrP}}JM zI9xU!a0XBd2m<5*1^^s@5CHM92@ne)`aJ+xz)J52krLl;0ua<^KLb_(h60FpqC@Rd zA6Egc19kw&q4xnC2Fw6F1t9)M0yYA^3p|sYiUZH}NU1K-p*rgUl>m}w5+DQ66W|5B z1n>c<{|dks0PW&`0Sp9m1W^5{02=({pQ$dT_>o_UpMN1GU3?87Igsp%05qwM1bfFvL&9L13U;I zhjTAL21o@EAGZNSm+I{WhyddIV}J&*4M1|EK9Szu23!GSpe}K z_=C@q?!>bG*-OVD&%8zOjzVg<=bL32^uDsYNTSnedYeYjd%}t<)HT+lZr)HHDal3- zAU@FISV#o;mxeR8D-D(dbDWsM*a@&}O6!XXw8l-C!Y*ym z+~~Xihtt;|Y4BY>{#a0_&sY6$!0jvNP|ndm`ip|NiGlWik$;lc$;6_*ofce%`C}>I gFMw5mX8~IPuLIr%d)3Hoe>iDR*RDbWk~PG$ z-<^?UIi!0yXLZigyt(hbcYoh~_uV=2HaBv>%0kTXH}%XKzPo3Y>MLN@_wn*++nyR` zeF^Pae1_k(uut)>eqy1_q}#4TBr}%#B*OF~Smt?cKNmmKHzu*Uv$_68>EfE%2qs!) z;ct{W)XaG9GlNC2+&;aBUS%iprQ;9gc{aqYLU`xk{CJWDu1uZ{T$(INabt6(piIy~6RMiv5nCGk+(cpEfyu*xtCJT4 z6_ZESRl1AUWS1WNv##EOHVCG+m_+Xj68O(e%`G75Frh(5+ zu)#d3T4AC#nqgTrcphZM;QACvR_#9+=Y+A7OteJx#8-7*;u`FG*4^t)-ffdQG4Qf(m5wVBc|g$x1NKR2;9 z>!G%HSHnKf7}mWdZCx&Kd(6aZ%)wik#aofZTam_Fo`bg}i`S6FTbRaks5ywWSwwpl zk&nGMQDvy=!>ce0imH)H58 zKG4ospiQ0P7i?u)-`bc_<@MFbjjX{2B~_9~&sO5?QsBZ^k#u~I4Pl!X!n@3{ZIx7ph4QflkM-n9AacuS&&))M<;=nm&)Eqg+}UHH zXqdA!!f1M2!1%0@EWE3?(mfNx4SYjbrH~t`9PY}X#C;HASYKQ8KHwvExQh7&MzE>Jf z9&OwP6i#>(y^Ne&47dQtF(ufHHhGp%^)Q~J28-Hm%5Zq)Zft|yU};X zTRPf^-Hr~q74LVZ$L_$K+f@w6n6yS(?|Ccv6~}r!$&K{HpQpE6`As8>zdy~Khm*@1 znbXz7oSK(8BVHE&gD$wqBZw(Ftoi6on^^pWP7w3>=A$aq8BCUeS zfJYDkNHcRmRMhEzxjvw<_}eMLUFl}=dff(Cot^|Nk`hz{SfS?wSgwEOA{!*4BOv$( zv^ZK4Z6>W6Rf5WKMLuj_NV;ETDc*Y&!NU#KCa0H^e$zc?^x@^7;|E> zPm@h_4bo)LHKj17A&_cdXiQ0g(QVyO;L|4uvyRPUagk~v|j3GkRm{;2BB|)iPODR?g535*J|dtIGXB> zV9v7fNbhMSTB#{j3PFbLb~hilX_{hfN-2RJ7szW~S1);uy$w&9`oZ2VFp@v!QqwQz z>@>s%wccvTNK;RHgWgClG~*+s+1v1xsUPu~Q!{ z7NLF=#uWJrq_>(_Qc?B9@H9CBDX1KzbN!&0)D6oAl-Cr>`R>vEUMvVm&{yfMOcSGg z*DsSVX#OfU`OKy10hx$#qs&L0mF7O+;QbJB)5!`!Ba|^+NQ}lj%8%>=3)k) zhQHPr{u;_*bGEyZrx2SN^8a=nmo?7|(=RlT?ye(&J~OQf?%&y6JuY&Q3(V<`1rT(| z6I{JPFh(vi2F`;$fvi^q(_V27UJ>Yh4ZDJmm#55e+bD#Abx$ipl>t0WZ^%n&N-<2D z@K6hBJ(6%EDOCSIvRNfc6-9 zIx#w%%{L7#U60i?o2opctlAu9RYO^8#7Zcs{}Rg~>hj^R@{RpPy^LR1FS$A~4O5q< zrIKZ8KDon=Rjk>yt}|v1X`WlgY2fjEGq6Zxq)K zMPSBizaOOx^i*`Wi*o_wo0?%uJyDP7MQ?PzwGH8sYqMaqzHsTy=lWH( z+Sv1-RqgFoTiZQzOWQdtbj*&j(Tah|2T;YCYg{Rxo0x+u;n_g32S*2X#9I?A5QJ>B z4f{U~s3T*P^P9JYZ!OqUI`R6y3!!gUN>qFp00T3G7lH8iP(vw`UeNfmvEuolC|p*N zKA}=&lAQOcNRH)@WP=BPUA-Z6sE*O)bT+{n5PT2>OAUfuMXeu{>?4x76pK`YBwZ=5 zC*7}eNPVv2$3s*km`#K(+Ur1M{wIiVO*4WO8&ry}Lq#Pej^~hJgR=>lr#o%ZqZ_Hz z5Rp+P!-ZLsOJnmqn!)xkjN(%0(tdO!-RI9v@PPDc?hNi>|`fq_bpJ&0^3jSg;c|C_xVCpipNf7?tfgAIqPC1due+q^zo&XA8u?Wdh=s*4Mx>s$p)p+ z*65IJQcHs!B^OFZP)-M@GkHrM80an~O_9d<7GlOswEe~rMyWNzt&wis603%eg2(_| z-YUuS(LZCQbxNO70VWn}NlJLc1!MG|rw-|f@1opF=yPK0uxW?Gr#mBSWY7wk(A=_qg^T~@WB89erO*-dlHp;Q3mkn5 z@Izn~%YlRy(VODAj!b~6@h^NB+HcmYX0n^^`OIOa(L11K*FLYICtJM?scyvWmJ{wS(wCBAEQn9`T;`1(C9sWA?OY8+Y|lG0nA4 zPXwo1fc5>9tIrpU;ArPAxcYXY7wlO}TI?-r>;F&-ZNXGey7?=FZ1 z6;B0E8t@&clG98@r~!Xs|B7iKUuO2WUMGWW}^W=e2LU^}C3ypO+tYeX$H`bYPX+%(3jP^Dy7=F3FM$i*=;h)!w z7A~%(!m0WOi=H??6?Qb(!(VII+Qf^o^HyAD%zgBGh#XzTtdFf`*5zn7pk?-7v)pX# zyGO%C&SsS{^9D!3qTrHMsB*7W^u#Zbd{G$ShIr?_G0N3Bq}yV95_|ViY$>Yj>4}du zzjtM9sR4a33)KN_h&32cZx+gb?Dtm0Y7OXTve3K=fQw^Q23X31tsjSrmsW=Bmdu6J zT==FXRj9l+LWWgnYtdTJHldN6oAD(1t|v#w{~?Dqm678=bL5C;p+C)$=d_F*FxKPAaO09n$6pp?n;}T}yUW_R_2G-l_Neu^jxCHve2eg0 zIuh~K;fdtjSBqzPv^GZennJ8ariB{B;iR*%G6!}f4V#+-;}P}<&k`$yXBTfex{(b& zgZBb4Z-kzOo_p|*s18qBt{Qa|Qt`k3aNzJuer3)Jcy@!!M)>9Y&1HAx#a6f%)H>ZY zweuCbT!KRHN5k_r^O)|GEvx3@eQJ0P4*9)2Rck3BHoqhcUVaFlfu|EejSXH!X^~;Q zSDTrYKAKYd18A~|O7kNHp;%9?oZz3t=c=!<_!m;)|6WnIkc)pO+2F(%V4EPH`5*pC z>lrTotz_6H+WDtfhnI>sVxxR9*}q+Up36T9o?l!Yp1xsK+pQ@U7Z|tvA10V>9N(Nv zqWtp0=#kf0d{t!(ITyY*?Vzq|QoJpZ#~ zc=`I<`JGq$mGv!L`Py5V_0Uwq9RuK=z3}*rx7_^H;(o`*O0JB>2U1t?Nt_hom|O;S z>-wclU7Wz;T==E-{rqFi;W|$zUXO-@o`ZPa+s{cq;4NFBkF_SEQTyJF;a9fa5Pos{_6lX!gL`2w zwqJSZQO5zrx8Gqv!Udh{!ax0-BfPA$BD}XVTzPYI%Poz!+`6)H?aIbgHcW7=bTB-_ zQrD8aZS8mc)otOYx||C?MitqDNL7wjgSHuM+p6%dyO#Dpecv%5{AI6Ot^KWj;kDg! z6BYO0{c`k)2Ry(3?MH_OrgW9VstbP6RluPmmqXD!F*fs)anf81JB0PWf_4P$DBAOA ZFQW~hy@AH&2iSHnn^{Ex$6I>t{9iYW=*j>9 diff --git a/G4:G431/CORDIC/cordic.c b/G4:G431/CORDIC/cordic.c index 1d74445..359de15 100644 --- a/G4:G431/CORDIC/cordic.c +++ b/G4:G431/CORDIC/cordic.c @@ -33,17 +33,12 @@ static void cordic_enable(void){ static int32_t float_to_q31(float x){ if(x >= 1.0f) x = 1.0f - 1e-6f; if(x <= -1.0f) x = -1.0f + 1e-6f; - return (int32_t)(x * 2147483648.0f); + return (int32_t)(x * Q31_BASE); } // Convert Q1.31 to float static float q31_to_float(int32_t q){ - return (float)q / 2147483648.0f; -} - -// Wait CORDIC ready -static void cordic_wait_ready(void){ - while(!(CORDIC->CSR & CORDIC_CSR_RRDY)); + return (float)q / Q31_BASE; } static void cordic_init(void){ @@ -72,28 +67,40 @@ static float cordic_scalar(uint32_t func_mode, float x, int arg_count){ } #endif +// sincos +void cordic_sincos(float angle, float *s, float *c){ + float norm = angle / M_PIf; + if(norm > 1.0f) norm = 1.0f; + if(norm < -1.0f) norm = -1.0f; + cordic_init(); + CORDIC->CSR = (CORDIC_CSR_FUNC_SIN << CORDIC_CSR_FUNC_Pos) | CORDIC_CSR_DEF | CORDIC_CSR_NRES; + CORDIC->WDATA = float_to_q31(norm); + int32_t res = CORDIC->RDATA; + if(s) *s = q31_to_float(res); + res = CORDIC->RDATA; + if(c) *c = q31_to_float(res); +} + // sin float cordic_sin(float angle){ - float norm = angle / 3.141592653589793f; + float norm = angle / M_PIf; if(norm > 1.0f) norm = 1.0f; if(norm < -1.0f) norm = -1.0f; cordic_init(); CORDIC->CSR = (CORDIC_CSR_FUNC_SIN << CORDIC_CSR_FUNC_Pos) | CORDIC_CSR_DEF; CORDIC->WDATA = float_to_q31(norm); - cordic_wait_ready(); - int32_t res = CORDIC->RDATA; // ÐÅÒ×ÏÅ ÞÔÅÎÉÅ -> sin + int32_t res = CORDIC->RDATA; return q31_to_float(res); } // cos float cordic_cos(float angle){ - float norm = angle / 3.141592653589793f; + float norm = angle / M_PIf; if(norm > 1.0f) norm = 1.0f; if(norm < -1.0f) norm = -1.0f; cordic_init(); CORDIC->CSR = (CORDIC_CSR_FUNC_COS << CORDIC_CSR_FUNC_Pos) | CORDIC_CSR_DEF; CORDIC->WDATA = float_to_q31(norm); - cordic_wait_ready(); int32_t res = CORDIC->RDATA; // cos return q31_to_float(res); } @@ -104,9 +111,8 @@ float cordic_atan(float val){ cordic_init(); CORDIC->CSR = (CORDIC_CSR_FUNC_ATAN << CORDIC_CSR_FUNC_Pos) | CORDIC_CSR_DEF; CORDIC->WDATA = float_to_q31(val); - cordic_wait_ready(); int32_t res = CORDIC->RDATA; - return q31_to_float(res) * 3.141592653589793f; + return q31_to_float(res) * M_PIf; } // sqrt @@ -123,7 +129,6 @@ float cordic_sqrt(float x){ cordic_init(); CORDIC->CSR = (CORDIC_CSR_FUNC_SQRT << CORDIC_CSR_FUNC_Pos) | CORDIC_CSR_DEF; CORDIC->WDATA = float_to_q31(x); - cordic_wait_ready(); int32_t res = CORDIC->RDATA; return scale * q31_to_float(res); } @@ -141,7 +146,6 @@ float cordic_log(float x){ cordic_init(); CORDIC->CSR = (CORDIC_CSR_FUNC_LOG << CORDIC_CSR_FUNC_Pos) | CORDIC_CSR_DEF; CORDIC->WDATA = float_to_q31(x); - cordic_wait_ready(); int32_t res = CORDIC->RDATA; return q31_to_float(res) * 0.6931471805599453f + add; } diff --git a/G4:G431/CORDIC/cordic.files b/G4:G431/CORDIC/cordic.files index a700edc..3f705b7 100644 --- a/G4:G431/CORDIC/cordic.files +++ b/G4:G431/CORDIC/cordic.files @@ -1,3 +1,5 @@ +astro.c +astro.h commproto.cpp commproto.h cordic.c diff --git a/G4:G431/CORDIC/cordic.h b/G4:G431/CORDIC/cordic.h index 05849bc..4d4e26d 100644 --- a/G4:G431/CORDIC/cordic.h +++ b/G4:G431/CORDIC/cordic.h @@ -20,6 +20,12 @@ #include +#ifndef M_PIf +#define M_PIf 3.141592653589793f +#endif + +#define Q31_BASE 2147483648.0f + // functions enum { CORDIC_CSR_FUNC_COS = 0, @@ -33,6 +39,7 @@ enum { CORDIC_CSR_FUNC_SQRT }; +void cordic_sincos(float angle, float *s, float *c); float cordic_sin(float angle); float cordic_cos(float angle); float cordic_atan(float val); diff --git a/G4:G431/CORDIC/strfunc.c b/G4:G431/CORDIC/strfunc.c index cc39bb7..809ac23 100644 --- a/G4:G431/CORDIC/strfunc.c +++ b/G4:G431/CORDIC/strfunc.c @@ -15,6 +15,7 @@ * along with this program. If not, see . */ +#include // isnan/isinf #include #include @@ -262,3 +263,116 @@ char *getint(char *txt, int32_t *I){ *I = sign * (int32_t)U; return nxt; } + +// be careful: if pow10 would be bigger you should change str[] size! +static const float pwr10[] = {1.f, 10.f, 100.f, 1000.f, 10000.f, 100000.f, 1000000.f, 10000000.f}; +static const float rounds[] = {0.5f, 0.05f, 0.005f, 0.0005f, 0.00005f, 0.000005f, 0.0000005f, 0.00000005f}; +#define P10L (sizeof(pwr10)/sizeof(float) - 1) +char *float2str(float x, uint8_t prec){ + static char str[16] = {0}; // -117.5494E-36\0 - 14 symbols max! + if(prec > P10L) prec = P10L; + if(isnan(x)){ memcpy(str, "NAN", 4); return str;} + else{ + int i = isinf(x); + if(i){memcpy(str, "-INF", 5); if(i == 1) return str+1; else return str;} + } + char *s = str + 14; // go to end of buffer + uint8_t minus = 0; + if(x < 0){ + x = -x; + minus = 1; + } + int pow = 0; // xxxEpow + // now convert float to 1.xxxE3y + while(x > 1000.f){ + x /= 1000.f; + pow += 3; + } + if(x > 0.) while(x < 1.){ + x *= 1000.f; + pow -= 3; + } + // print Eyy + if(pow){ + uint8_t m = 0; + if(pow < 0){pow = -pow; m = 1;} + while(pow){ + int p10 = pow/10; + *s-- = '0' + (pow - 10*p10); + pow = p10; + } + if(m) *s-- = '-'; + *s-- = 'E'; + } + // now our number is in [1, 1000] + uint32_t units; + if(prec){ + units = (uint32_t) x; + uint32_t decimals = (uint32_t)((x-units+rounds[prec])*pwr10[prec]); + // print decimals + while(prec){ + int d10 = decimals / 10; + *s-- = '0' + (decimals - 10*d10); + decimals = d10; + --prec; + } + // decimal point + *s-- = '.'; + }else{ // without decimal part + units = (uint32_t) (x + 0.5); + } + // print main units + if(units == 0) *s-- = '0'; + else while(units){ + uint32_t u10 = units / 10; + *s-- = '0' + (units - 10*u10); + units = u10; + } + if(minus) *s-- = '-'; + return s+1; +} + +char *getfloat(char *str, float *f){ + str = omit_spaces(str); + int isminus = 0; + if(*str == '-'){ + isminus = 1; + ++str; + }else if (*str == '+'){ + ++str; + } + float result = 0.0f; + while(*str >= '0' && *str <= '9'){ + result = result * 10.0f + (float)(*str - '0'); + ++str; + } + if(*str != '.') goto retn; + + ++str; // omit point + float frac = 0.0f; + float divisor = 1.0f; + + while(*str >= '0' && *str <= '9'){ + divisor *= 10.0f; + frac = frac * 10.0f + (float)(*str - '0'); + ++str; + } + + if(divisor > 1.0f){ + result += frac / divisor; + } + + if(*str == 'e' || *str == 'E'){ // exp + ++str; + int32_t E; + if(str != getint(str, &E)){ + if(E > 0) while(E-- > 0) result *= 10.f; + else while(E++ < 0) result /= 10.f; + } + } + +retn: + if(f) *f = (isminus) ? -result : result; + return str; +} + diff --git a/G4:G431/CORDIC/strfunc.h b/G4:G431/CORDIC/strfunc.h index bd1950c..1e7cce2 100644 --- a/G4:G431/CORDIC/strfunc.h +++ b/G4:G431/CORDIC/strfunc.h @@ -19,12 +19,13 @@ #include +char *omit_spaces(char *buf); void hexdump(int (*sendfun)(const char*), uint8_t *arr, uint16_t len); char *u2str(uint32_t val); char *i2str(int32_t i); +char *float2str(float x, uint8_t prec); char *uhex2str(uint32_t val); char *gethex(const char *buf, uint32_t *N); char *getnum(char *txt, uint32_t *N); -char *omit_spaces(char *buf); char *getint(char *txt, int32_t *I); - +char *getfloat(char *str, float *f); diff --git a/G4:G431/CORDIC/test.c b/G4:G431/CORDIC/test.c index aae3602..12472f0 100644 --- a/G4:G431/CORDIC/test.c +++ b/G4:G431/CORDIC/test.c @@ -32,7 +32,7 @@ static float arr[N_TESTS]; // RNG static uint32_t rand_state = 123456789; -static uint32_t next_rand(void){ +static uint32_t next_rand(){ rand_state = rand_state * 1664525 + 1013904223; return rand_state; } @@ -76,36 +76,52 @@ static uint32_t run_test(void (*gen)(), float (*func)(float)){ return timer_read(); } +static uint32_t run_test2(void (*gen)(), void (*func)(float, float*, float*)){ + gen(); + volatile float result1 = 0.f, result2 = 0.f; // don't let gcc to optimize this cycle + timer_start(); + for(int i = 0; i < N_TESTS; ++i){ + func(arr[i], (float*)&result1, (float*)&result2); + (void) result1; + (void) result2; + } + timer_stop(); + return timer_read(); +} + // ------------- math.h tests ------------- -uint32_t test_math_sin(void){ +uint32_t test_math_sin(){ return run_test(fill_random_sin_cos, sinf); } -uint32_t test_math_cos(void){ +uint32_t test_math_cos(){ return run_test(fill_random_sin_cos, cosf); } -uint32_t test_math_atan(void){ +uint32_t test_math_atan(){ return run_test(fill_random_atan, atanf); } -uint32_t test_math_sqrt(void){ +uint32_t test_math_sqrt(){ return run_test(fill_random_sqrt, sqrtf); } -uint32_t test_math_log(void){ +uint32_t test_math_log(){ return run_test(fill_random_log, logf); } // ------------- CORDIC tests ------------- -uint32_t test_cordic_sin(void){ +uint32_t test_cordic_sincos(){ + return run_test2(fill_random_sin_cos, cordic_sincos); +} +uint32_t test_cordic_sin(){ return run_test(fill_random_sin_cos, cordic_sin); } -uint32_t test_cordic_cos(void){ +uint32_t test_cordic_cos(){ return run_test(fill_random_sin_cos, cordic_cos); } -uint32_t test_cordic_atan(void){ +uint32_t test_cordic_atan(){ return run_test(fill_random_atan, cordic_atan); } -uint32_t test_cordic_sqrt(void){ +uint32_t test_cordic_sqrt(){ return run_test(fill_random_sqrt, cordic_sqrt); } -uint32_t test_cordic_log(void){ +uint32_t test_cordic_log(){ return run_test(fill_random_log, cordic_log); } diff --git a/G4:G431/CORDIC/test.h b/G4:G431/CORDIC/test.h index fe6d095..b84248a 100644 --- a/G4:G431/CORDIC/test.h +++ b/G4:G431/CORDIC/test.h @@ -28,6 +28,7 @@ uint32_t test_math_sqrt(); uint32_t test_math_log(); // CORDIC tests +uint32_t test_cordic_sincos(); uint32_t test_cordic_sin(); uint32_t test_cordic_cos(); uint32_t test_cordic_atan(); diff --git a/G4:G431/CORDIC/version.inc b/G4:G431/CORDIC/version.inc index d8c9689..7630c46 100644 --- a/G4:G431/CORDIC/version.inc +++ b/G4:G431/CORDIC/version.inc @@ -1,2 +1,2 @@ -#define BUILD_NUMBER "19" -#define BUILD_DATE "2026-08-04" +#define BUILD_NUMBER "48" +#define BUILD_DATE "2026-08-12"