mirror
/
RS-tracker
kopia lustrzana https://github.com/rs1729/RS
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

rs41: RH2.1

pull/37/head
Zilog80 2021-01-23 20:01:02 +01:00
rodzic 3d33874049
commit 107360c68a
1 zmienionych plików z 112 dodań i 25 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

137
demod/mod/rs41mod.c
Wyświetl plik

@ -113,6 +113,8 @@ typedef struct {
float ptu_calT2[3]; // calibration T2-Hum
float ptu_calH[2]; // calibration Hum
float ptu_mtxH[42];
float ptu_corHp[3];
float ptu_corHt[12];
float ptu_Cf1;
float ptu_Cf2;
float ptu_calP[25];
@ -521,6 +523,12 @@ static int get_CalData(gpx_t *gpx) {
for (j = 0; j < 42; j++) { // 0x07*0x10+13 = 0x07D = 125
memcpy(gpx->ptu_mtxH+j, gpx->calibytes+125+4*j, 4);
}
for (j = 0; j < 3; j++) { // 0x2A*0x10+6 = 0x2A6 = 678
memcpy(gpx->ptu_corHp+j, gpx->calibytes+678+4*j, 4);
}
for (j = 0; j < 12; j++) { // 0x2B*0x10+A = 0x2BA = 698
memcpy(gpx->ptu_corHt+j, gpx->calibytes+698+4*j, 4);
}
// cf. DF9DQ or stsst/RS-fork
memcpy(gpx->ptu_calP+0, gpx->calibytes+606, 4); // 0x25*0x10+14 = 0x25E
memcpy(gpx->ptu_calP+4, gpx->calibytes+610, 4); // ..
@ -562,7 +570,7 @@ static float get_T(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float *ptu_co, fl
// (data:) ftp://ftp-cdc.dwd.de/climate_environment/CDC/observations_germany/radiosondes/
// (diffAlt: Ellipsoid-Geoid)
// (note: humidity sensor has significant time-lag at low temperatures)
static float get_RH(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T) {
static float get_RHemp(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T) {
float a0 = 7.5; // empirical
float a1 = 350.0/gpx->ptu_calH[0]; // empirical
float fh = (f-f1)/(float)(f2-f1);
@ -598,8 +606,8 @@ static float vaporSatP(float Tc) {
return (float)p; // [Pa]
}
// cf. github DF9DQ // offset stratosphere RH2(-60C)=+5% , RH2(-40C)=+1% ?
static float get_RH2adv(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T, float TH) {
// cf. github DF9DQ
static float get_RH2adv(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T, float TH, float P) {
float rh = 0.0;
float cfh = (f-f1)/(float)(f2-f1);
float cap = gpx->ptu_Cf1+(gpx->ptu_Cf2-gpx->ptu_Cf1)*cfh;
@ -612,10 +620,33 @@ static float get_RH2adv(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T, flo
bk = 1.0;
for (k = 0; k < 6; k++) {
b[k] = bk;
b[k] = bk; // Tp^k
bk *= Trh_20_180;
}
if (P > 0.0) // in particular if P<200hPa , T<-40
{
double _p = P / 1000.0; // bar
double _cpj = 1.0;
double corrCp = 0.0;
double bt, bp[3];
for (j = 0; j < 3; j++) {
bp[j] = gpx->ptu_corHp[j] * (_p/(1.0 + gpx->ptu_corHp[j]*_p) - _cpj/(1.0 + gpx->ptu_corHp[j]));
_cpj *= Cp; // Cp^j
}
corrCp = 0.0;
for (j = 0; j < 3; j++) {
bt = 0.0;
for (k = 0; k < 4; k++) {
bt += gpx->ptu_corHt[4*j+k] * b[k];
}
corrCp += bp[j] * bt;
}
Cp -= corrCp;
}
aj = 1.0;
for (j = 0; j < 7; j++) {
for (k = 0; k < 6; k++) {
@ -624,7 +655,7 @@ static float get_RH2adv(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T, flo
aj *= Cp;
}
if ( 1 ) { // empirical correction
if ( P <= 0.0 ) { // empirical correction
float T2 = -40;
if (T < T2) { _rh += (T-T2)/12.0; }
}
@ -659,9 +690,48 @@ static float get_P(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, int fx)
return (float)p;
}
// ---------------------------------------------------------------------------------------
//
// barometric formula https://en.wikipedia.org/wiki/Barometric_formula
static float Ph(float h) {
double Pb, Tb, Lb, hb;
//double RgM = 8.31446/(9.80665*0.0289644);
double gMR = 9.80665*0.0289644/8.31446;
float P = 0.0;
if (h > 32000.0) { //P < 8.6802
Pb = 8.6802;
Tb = 228.65;
Lb = 0.0028;
hb = 32000.0;
}
else if (h > 20000.0) { // P < 54.7489 (&& P >= 8.6802)
Pb = 54.7489;
Tb = 216.65;
Lb = 0.001;
hb = 20000.0;
}
else if (h > 11000.0) { // P < 226.321 (&& P >= 54.7489)
Pb = 226.321;
Tb = 216.65;
Lb = 0.0;
hb = 11000.0;
}
else { // P >= 226.321
Pb = 1013.25;
Tb = 288.15;
Lb = -0.0065;
hb = 0.0;
}
static int get_PTU(gpx_t *gpx, int ofs, int pck) {
//if (Lb == 0.0) altP = -RgM*Tb * log(P/Pb) + hb;
//else altP = Tb/Lb * (pow(P/Pb, -RgM*Lb)-1.0) + hb;
if (Lb == 0.0) P = Pb * exp( -gMR*(h-hb)/Tb );
else P = Pb * pow( 1.0+Lb*(h-hb)/Tb , -gMR/Lb);
return P;
}
static int get_PTU(gpx_t *gpx, int ofs, int pck, int valid_alt) {
int err=0, i;
int bR, bc1, bT1,
bc2, bT2;
@ -697,7 +767,10 @@ static int get_PTU(gpx_t *gpx, int ofs, int pck) {
bH2 = gpx->calfrchk[0x08] && gpx->calfrchk[0x09]
&& gpx->calfrchk[0x10] && gpx->calfrchk[0x11]
&& gpx->calfrchk[0x12] && gpx->calfrchk[0x13];
&& gpx->calfrchk[0x12] && gpx->calfrchk[0x13]
&& gpx->calfrchk[0x2A] && gpx->calfrchk[0x2B]
&& gpx->calfrchk[0x2C] && gpx->calfrchk[0x2D]
&& gpx->calfrchk[0x2E];
bP = gpx->calfrchk[0x21] && gpx->calibytes[0x21F] == 'P'
&& gpx->calfrchk[0x25] && gpx->calfrchk[0x26]
@ -715,21 +788,26 @@ static int get_PTU(gpx_t *gpx, int ofs, int pck) {
gpx->TH = TH;
if (bH && Tc > -273.0) {
RH = get_RH(gpx, meas[3], meas[4], meas[5], Tc); // TH, TH-Tc (sensorT - T)
RH = get_RHemp(gpx, meas[3], meas[4], meas[5], Tc); // TH, TH-Tc (sensorT - T)
}
gpx->RH = RH;
// cf. DF9DQ, stsst/RS-fork
if (gpx->option.ptu == 2) {
if (bH && bH2 && Tc > -273.0 && TH > -273.0) {
RH2 = get_RH2adv(gpx, meas[3], meas[4], meas[5], Tc, TH);
}
}
gpx->RH2 = RH2;
// cf. DF9DQ, stsst/RS-fork
if (bP) {
P = get_P(gpx, meas[9], meas[10], meas[11], i2(gpx->frame+pos_PTU+ofs+2+38));
}
gpx->P = P;
if (gpx->option.ptu == 2) {
float _P = -1.0;
if (bP) _P = P;
else { // approx
if (valid_alt > 0) _P = Ph(gpx->alt);
}
if (bH && bH2 && Tc > -273.0 && TH > -273.0) {
RH2 = get_RH2adv(gpx, meas[3], meas[4], meas[5], Tc, TH, _P);
}
}
gpx->RH2 = RH2;
if (gpx->option.vbs == 4 && (gpx->crc & (crc_PTU | crc_GPS3))==0)
@ -1447,7 +1525,7 @@ static int prn_frm(gpx_t *gpx) {
static int prn_ptu(gpx_t *gpx) {
fprintf(stdout, " ");
if (gpx->T > -273.0) fprintf(stdout, " T=%.1fC ", gpx->T);
if (gpx->RH > -0.5 && gpx->option.ptu != 2) fprintf(stdout, " RH=%.0f%% ", gpx->RH);
if (gpx->RH > -0.5 && gpx->option.ptu != 2) fprintf(stdout, " _RH=%.0f%% ", gpx->RH);
if (gpx->P > 0.0) {
if (gpx->P < 100.0) fprintf(stdout, " P=%.2fhPa ", gpx->P);
else fprintf(stdout, " P=%.1fhPa ", gpx->P);
@ -1574,6 +1652,7 @@ static int print_position(gpx_t *gpx, int ec) {
int out = 1;
int sat = 0;
int pos_aux = 0, cnt_aux = 0;
int ofs_ptu = 0, pck_ptu = 0;
//gpx->out = 0;
gpx->aux = 0;
@ -1627,10 +1706,11 @@ static int print_position(gpx_t *gpx, int ec) {
break;
case pck_PTU: // 0x7A2A
ofs = pos - pos_PTU;
err0 = get_PTU(gpx, ofs, pck_PTU);
if ( 0 && !err0 && gpx->option.ptu ) {
prn_ptu(gpx);
ofs_ptu = pos - pos_PTU;
pck_ptu = pck_PTU;
if ( 0 && gpx->option.ptu ) {
//err0 = get_PTU(gpx, ofs_ptu, pck_ptu);
// if (!err0) prn_ptu(gpx);
}
break;
@ -1662,8 +1742,11 @@ static int print_position(gpx_t *gpx, int ec) {
break;
case pck_SGM_xTU: // 0x7F1B
ofs = pos - pos_PTU;
err0 = get_PTU(gpx, ofs, pck);
ofs_ptu = pos - pos_PTU;
pck_ptu = pck;
if ( 0 ) {
//err0 = get_PTU(gpx, ofs_ptu, pck_ptu);
}
break;
case pck_SGM_CRYPT: // 0x80A7
@ -1695,9 +1778,11 @@ static int print_position(gpx_t *gpx, int ec) {
if ( pos > frm_end ) // end of (sub)frame
{
if (gpx->option.ptu && out && !sat && !err0 && !encrypted) {
prn_ptu(gpx);
if (gpx->option.ptu && out && !sat && !encrypted && pck_ptu > 0) {
err0 = get_PTU(gpx, ofs_ptu, pck_ptu, !err3);
if (!err0) prn_ptu(gpx);
}
pck_ptu = 0;
get_Calconf(gpx, out, ofs_cal);
@ -1774,7 +1859,7 @@ static int print_position(gpx_t *gpx, int ec) {
ofs = 0;
if (pck < 0x8000) {
err0 = get_PTU(gpx, 0, pck);
//err0 = get_PTU(gpx, 0, pck, 0);
if (pck == pck_PTU) ofs = 0;
else if (pck == pck_SGM_xTU) ofs = 0x1B-0x2A;
@ -1783,6 +1868,8 @@ static int print_position(gpx_t *gpx, int ec) {
err3 = get_GPS3(gpx, ofs);
if (!err1) Gps2Date(gpx);
err0 = get_PTU(gpx, 0, pck, !err3);
if (out) {
if (!err1) prn_gpstime(gpx);