multi: rs41 PTU: Pressure adopted from DF9DQ and stsst/RS-fork, more detailed rel.hum. cf. DF9DQ

pull/52/head
Zilog80 2020-12-22 19:09:35 +01:00
rodzic 6ba3386662
commit 7f3f0f4106
1 zmienionych plików z 201 dodań i 62 usunięć

Wyświetl plik

@ -77,7 +77,8 @@ typedef struct {
int std; int min; float sek;
double lat; double lon; double alt;
double vH; double vD; double vV;
float T; float RH;
float T; float RH; float TH;
float P; float RH2;
ui32_t crc;
ui8_t frame[FRAME_LEN];
ui8_t calibytes[51*16];
@ -89,6 +90,10 @@ typedef struct {
float ptu_co2[3]; // { -243.911 , 0.187654 , 8.2e-06 }
float ptu_calT2[3]; // calibration T2-Hum
float ptu_calH[2]; // calibration Hum
float ptu_mtxH[42];
float ptu_Cf1;
float ptu_Cf2;
float ptu_calP[25];
ui32_t freq; // freq/kHz (RS41)
int jsn_freq; // freq/kHz (SDR)
float batt; // battery voltage (V)
@ -203,6 +208,13 @@ static ui32_t u2(ui8_t *bytes) { // 16bit unsigned int
return bytes[0] | (bytes[1]<<8);
}
static int i2(ui8_t *bytes) { // 16bit signed int
//return (i16_t)u2(bytes);
int val = bytes[0] | (bytes[1]<<8);
if (val & 0x8000) val -= 0x10000;
return val;
}
/*
double r8(ui8_t *bytes) {
double val = 0;
@ -363,14 +375,14 @@ static int get_BattVolts(gpx_t *gpx, int ofs) {
int i;
unsigned byte;
ui8_t batt_bytes[2];
float batt_volts;
ui16_t batt_volts; // signed voltage?
for (i = 0; i < 2; i++) {
byte = gpx->frame[pos_BattVolts+ofs + i];
batt_bytes[i] = byte;
}
batt_volts = (float)(batt_bytes[0] + (batt_bytes[1] << 8));
// 2 bytes? V > 25.5 ?
batt_volts = batt_bytes[0]; // + (batt_bytes[1] << 8);
gpx->batt = batt_volts/10.0;
return 0;
@ -402,6 +414,8 @@ static int get_SondeID(gpx_t *gpx, int crc, int ofs) {
// don't reset gpx->frame[] !
gpx->T = -273.15;
gpx->RH = -1.0;
gpx->P = -1.0;
gpx->RH2 = -1.0;
// new ID:
memcpy(gpx->id, sondeid_bytes, 8);
gpx->id[8] = '\0';
@ -440,6 +454,8 @@ static int get_FrameConf(gpx_t *gpx, int ofs) {
static int get_CalData(gpx_t *gpx) {
int j;
memcpy(&(gpx->ptu_Rf1), gpx->calibytes+61, 4); // 0x03*0x10+13
memcpy(&(gpx->ptu_Rf2), gpx->calibytes+65, 4); // 0x04*0x10+ 1
@ -450,6 +466,7 @@ static int get_CalData(gpx_t *gpx) {
memcpy(gpx->ptu_calT1+0, gpx->calibytes+89, 4); // 0x05*0x10+ 9
memcpy(gpx->ptu_calT1+1, gpx->calibytes+93, 4); // 0x05*0x10+13
memcpy(gpx->ptu_calT1+2, gpx->calibytes+97, 4); // 0x06*0x10+ 1
// ptu_calT1[3..6]
memcpy(gpx->ptu_calH+0, gpx->calibytes+117, 4); // 0x07*0x10+ 5
memcpy(gpx->ptu_calH+1, gpx->calibytes+121, 4); // 0x07*0x10+ 9
@ -461,82 +478,168 @@ static int get_CalData(gpx_t *gpx) {
memcpy(gpx->ptu_calT2+0, gpx->calibytes+305, 4); // 0x13*0x10+ 1
memcpy(gpx->ptu_calT2+1, gpx->calibytes+309, 4); // 0x13*0x10+ 5
memcpy(gpx->ptu_calT2+2, gpx->calibytes+313, 4); // 0x13*0x10+ 9
// ptu_calT2[3..6]
// cf. DF9DQ
memcpy(&(gpx->ptu_Cf1), gpx->calibytes+69, 4); // 0x04*0x10+ 5
memcpy(&(gpx->ptu_Cf2), gpx->calibytes+73, 4); // 0x04*0x10+ 9
for (j = 0; j < 42; j++) { // 0x07*0x10+13 = 0x07D = 125
memcpy(gpx->ptu_mtxH+j, gpx->calibytes+125+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); // ..
memcpy(gpx->ptu_calP+8, gpx->calibytes+614, 4);
memcpy(gpx->ptu_calP+12, gpx->calibytes+618, 4);
memcpy(gpx->ptu_calP+16, gpx->calibytes+622, 4);
memcpy(gpx->ptu_calP+20, gpx->calibytes+626, 4);
memcpy(gpx->ptu_calP+24, gpx->calibytes+630, 4);
memcpy(gpx->ptu_calP+1, gpx->calibytes+634, 4);
memcpy(gpx->ptu_calP+5, gpx->calibytes+638, 4);
memcpy(gpx->ptu_calP+9, gpx->calibytes+642, 4);
memcpy(gpx->ptu_calP+13, gpx->calibytes+646, 4);
memcpy(gpx->ptu_calP+2, gpx->calibytes+650, 4);
memcpy(gpx->ptu_calP+6, gpx->calibytes+654, 4);
memcpy(gpx->ptu_calP+10, gpx->calibytes+658, 4);
memcpy(gpx->ptu_calP+14, gpx->calibytes+662, 4);
memcpy(gpx->ptu_calP+3, gpx->calibytes+666, 4);
memcpy(gpx->ptu_calP+7, gpx->calibytes+670, 4); // ..
memcpy(gpx->ptu_calP+11, gpx->calibytes+674, 4); // 0x2A*0x10+ 2
return 0;
}
/*
static float get_Tc0(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2) {
// y = (f - f1) / (f2 - f1);
// y1 = (f - f1) / f2; // = (1 - f1/f2)*y
float a = 3.9083e-3, // Pt1000 platinum resistance
b = -5.775e-7,
c = -4.183e-12; // below 0C, else C=0
float *cal = gpx->ptu_calT1;
float Rb = (f1*gpx->ptu_Rf2-f2*gpx->ptu_Rf1)/(f2-f1), // ofs
Ra = f * (gpx->ptu_Rf2-gpx->ptu_Rf1)/(f2-f1) - Rb,
raw = Ra/1000.0,
g_r = 0.8024*cal[0] + 0.0176, // empirisch
r_o = 0.0705*cal[1] + 0.0011, // empirisch
r = raw * g_r + r_o,
t = (-a + sqrt(a*a + 4*b*(r-1)))/(2*b); // t>0: c=0
// R/R0 = 1 + at + bt^2 + c(t-100)t^3 , R0 = 1000 Ohm, t/Celsius
return t;
}
*/
// T_RH-sensor
static float get_TH(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2) {
float *p = gpx->ptu_co2;
float *c = gpx->ptu_calT2;
// temperature, platinum resistor
// T-sensor: gpx->ptu_co1 , gpx->ptu_calT1
// T_RH-sensor: gpx->ptu_co2 , gpx->ptu_calT2
static float get_T(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float *ptu_co, float *ptu_calT) {
float *p = ptu_co;
float *c = ptu_calT;
float g = (float)(f2-f1)/(gpx->ptu_Rf2-gpx->ptu_Rf1), // gain
Rb = (f1*gpx->ptu_Rf2-f2*gpx->ptu_Rf1)/(float)(f2-f1), // ofs
Rc = f/g - Rb,
//R = (Rc + c[1]) * c[0],
//T = p[0] + p[1]*R + p[2]*R*R;
R = Rc * c[0],
T = (p[0] + p[1]*R + p[2]*R*R + c[1])*(1.0 + c[2]);
return T;
}
// T-sensor, platinum resistor
static float get_Tc(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2) {
float *p = gpx->ptu_co1;
float *c = gpx->ptu_calT1;
float g = (float)(f2-f1)/(gpx->ptu_Rf2-gpx->ptu_Rf1), // gain
Rb = (f1*gpx->ptu_Rf2-f2*gpx->ptu_Rf1)/(float)(f2-f1), // ofs
Rc = f/g - Rb,
//R = (Rc + c[1]) * c[0],
//T = p[0] + p[1]*R + p[2]*R*R;
R = Rc * c[0],
T = (p[0] + p[1]*R + p[2]*R*R + c[1])*(1.0 + c[2]);
return T;
return T; // [Celsius]
}
// rel.hum., capacitor
// (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) {
float a0 = 7.5; // empirical
float a1 = 350.0/gpx->ptu_calH[0]; // empirical
float fh = (f-f1)/(float)(f2-f1);
float rh = 100.0 * ( a1*fh - a0 );
float T0 = 0.0, T1 = -25.0; // T/C
rh += T0 - T/5.5; // empir. temperature compensation
if (T < T1) rh *= 1.0 + (T1-T)/90.0; // empir. temperature compensation
float T0 = 0.0, T1 = -20.0, T2 = -40.0; // T/C v0.4
rh += T0 - T/5.5; // empir. temperature compensation
if (T < T1) rh *= 1.0 + (T1-T)/100.0; // empir. temperature compensation
if (T < T2) rh *= 1.0 + (T2-T)/120.0; // empir. temperature compensation
if (rh < 0.0) rh = 0.0;
if (rh > 100.0) rh = 100.0;
if (T < -273.0) rh = -1.0;
return rh;
}
// ---------------------------------------------------------------------------------------
//
// cf. github DF9DQ
// water vapor saturation pressure (Hyland and Wexler)
static float vaporSatP(float Tc) {
double T = Tc + 273.15;
// Apply some correction magic
// T = -0.4931358 + (1.0 + 4.61e-3) * T - 1.3746454e-5 * T*T + 1.2743214e-8 * T*T*T;
// H+W equation
double p = expf(-5800.2206 / T
+1.3914993
+6.5459673 * log(T)
-4.8640239e-2 * T
+4.1764768e-5 * T*T
-1.4452093e-8 * T*T*T
);
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) {
float rh = 0.0;
float cfh = (f-f1)/(float)(f2-f1);
float cap = gpx->ptu_Cf1+(gpx->ptu_Cf2-gpx->ptu_Cf1)*cfh;
double Cp = (cap / gpx->ptu_calH[0] - 1.0) * gpx->ptu_calH[1];
double Trh_20_180 = (TH - 20.0) / 180.0;
double _rh = 0.0;
double aj = 1.0;
double bk = 1.0, b[6];
int j, k;
bk = 1.0;
for (k = 0; k < 6; k++) {
b[k] = bk;
bk *= Trh_20_180;
}
aj = 1.0;
for (j = 0; j < 7; j++) {
for (k = 0; k < 6; k++) {
_rh += aj * b[k] * gpx->ptu_mtxH[6*j+k];
}
aj *= Cp;
}
if ( 1 ) { // empirical correction
float T2 = -40;
if (T < T2) { _rh += (T-T2)/12.0; }
}
rh = _rh * vaporSatP(TH)/vaporSatP(T);
if (rh < 0.0) rh = 0.0;
if (rh > 100.0) rh = 100.0;
return rh;
}
//
// cf. github DF9DQ or stsst/RS-fork
static float get_P(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, int fx)
{
double p = 0.0;
double a0, a1, a0j, a1k;
int j, k;
if (f1 == f2 || f1 == f) return 0.0;
a0 = gpx->ptu_calP[24] / ((float)(f - f1) / (float)(f2 - f1));
a1 = fx * 0.01;
a0j = 1.0;
for (j = 0; j < 6; j++) {
a1k = 1.0;
for (k = 0; k < 4; k++) {
p += a0j * a1k * gpx->ptu_calP[j*4+k];
a1k *= a1;
}
a0j *= a0;
}
return (float)p;
}
// ---------------------------------------------------------------------------------------
static int get_PTU(gpx_t *gpx, int ofs, int pck) {
int err=0, i;
int bR, bc1, bT1,
bc2, bT2;
int bH;
int bH2;
int bP;
ui32_t meas[12];
float Tc = -273.15;
float TH = -273.15;
float RH = -1.0;
float RH2 = -1.0;
float P = -1.0;
get_CalData(gpx);
@ -558,21 +661,42 @@ static int get_PTU(gpx_t *gpx, int ofs, int pck) {
bT2 = gpx->calfrchk[0x13];
bH = gpx->calfrchk[0x07];
bH2 = gpx->calfrchk[0x08] && gpx->calfrchk[0x09]
&& gpx->calfrchk[0x10] && gpx->calfrchk[0x11]
&& gpx->calfrchk[0x12] && gpx->calfrchk[0x13];
bP = gpx->calfrchk[0x21] && gpx->calibytes[0x21F] == 'P'
&& gpx->calfrchk[0x25] && gpx->calfrchk[0x26]
&& gpx->calfrchk[0x27] && gpx->calfrchk[0x28]
&& gpx->calfrchk[0x29] && gpx->calfrchk[0x2A];
if (bR && bc1 && bT1) {
Tc = get_Tc(gpx, meas[0], meas[1], meas[2]);
//Tc0 = get_Tc0(gpx, meas[0], meas[1], meas[2]);
Tc = get_T(gpx, meas[0], meas[1], meas[2], gpx->ptu_co1, gpx->ptu_calT1);
}
gpx->T = Tc;
if (bR && bc2 && bT2) {
TH = get_TH(gpx, meas[6], meas[7], meas[8]);
TH = get_T(gpx, meas[6], meas[7], meas[8], gpx->ptu_co2, gpx->ptu_calT2);
}
gpx->TH = TH;
if (bH) {
if (bH && Tc > -273.0) {
RH = get_RH(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;
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.vbs == 4 && (gpx->crc & (crc_PTU | crc_GPS3))==0)
{
@ -585,7 +709,7 @@ static int get_PTU(gpx_t *gpx, int ofs, int pck) {
printf("3: %8d %8d %8d", meas[6], meas[7], meas[8]);
printf(" # ");
//if (Tc > -273.0 && RH > -0.5)
if (0 && Tc > -273.0 && RH > -0.5)
{
printf(" ");
printf(" Tc:%.2f ", Tc);
@ -696,6 +820,7 @@ static int get_GPS2(gpx_t *gpx, int ofs) {
return err;
}
// WGS84/GRS80 Ellipsoid
#define EARTH_a 6378137.0
#define EARTH_b 6356752.31424518
#define EARTH_a2_b2 (EARTH_a*EARTH_a - EARTH_b*EARTH_b)
@ -889,10 +1014,10 @@ static int get_Calconf(gpx_t *gpx, int out, int ofs) {
byte = gpx->frame[pos_CalData+ofs+1+i];
fprintf(stdout, "%02x ", byte);
}
/*
/*
if (err == 0) fprintf(stdout, "[OK]");
else fprintf(stdout, "[NO]");
*/
*/
fprintf(stdout, " ");
}
@ -1074,6 +1199,13 @@ 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) 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);
}
if (gpx->option.ptu == 2) {
if (gpx->RH2 > -0.5) fprintf(stdout, " RH2=%.0f%% ", gpx->RH2);
}
return 0;
}
@ -1154,7 +1286,7 @@ static int prn_sat3(gpx_t *gpx, int ofs) {
pDOP = gpx->frame[pos_pDOP+ofs]/10.0; if (gpx->frame[pos_pDOP+ofs] == 0xFF) pDOP = -1.0;
fprintf(stdout, "numSatsFix: %2d sAcc: %.1f pDOP: %.1f\n", numSV, sAcc, pDOP);
/*
/*
fprintf(stdout, "CRC: ");
fprintf(stdout, " %04X", pck_GPS1);
if (check_CRC(gpx, pos_GPS1+ofs, pck_GPS1)==0) fprintf(stdout, "[OK]"); else fprintf(stdout, "[NO]");
@ -1167,12 +1299,12 @@ static int prn_sat3(gpx_t *gpx, int ofs) {
//fprintf(stdout, "[%+d]", check_CRC(gpx, pos_GPS3, pck_GPS3));
fprintf(stdout, "\n");
*/
*/
return 0;
}
static int print_position(gpx_t *gpx, int ec) {
int i, j;
int i;
int err, err0, err1, err2, err3;
//int output, out_mask;
int encrypted = 0;
@ -1325,11 +1457,18 @@ static int print_position(gpx_t *gpx, int ec) {
fprintf(stdout, "{ \"type\": \"%s\"", "RS41");
fprintf(stdout, ", \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f, \"sats\": %d, \"bt\": %d, \"batt\": %.2f",
gpx->frnr, gpx->id, gpx->jahr, gpx->monat, gpx->tag, gpx->std, gpx->min, gpx->sek, gpx->lat, gpx->lon, gpx->alt, gpx->vH, gpx->vD, gpx->vV, gpx->numSV, gpx->conf_cd, gpx->batt );
if (gpx->option.ptu && !err0 && gpx->T > -273.0) {
fprintf(stdout, ", \"temp\": %.1f", gpx->T );
}
if (gpx->option.ptu && !err0 && gpx->RH > -0.5) {
fprintf(stdout, ", \"humidity\": %.1f", gpx->RH );
if (gpx->option.ptu && !err0) {
float _RH = gpx->RH;
if (gpx->option.ptu == 2) _RH = gpx->RH2;
if (gpx->T > -273.0) {
fprintf(stdout, ", \"temp\": %.1f", gpx->T );
}
if (_RH > -0.5) {
fprintf(stdout, ", \"humidity\": %.1f", _RH );
}
if (gpx->P > 0.0) {
fprintf(stdout, ", \"pressure\": %.2f", gpx->P );
}
}
if (gpx->aux) { // <=> gpx->xdata[0]!='\0'
fprintf(stdout, ", \"aux\": \"%s\"", gpx->xdata );
@ -1523,7 +1662,7 @@ void *thd_rs41(void *targs) { // pcm_t *pcm, double xlt_fq
// init gpx
gpx.option.vbs = 1;
gpx.option.ptu = 1;
gpx.option.ptu = 2;
gpx.option.aut = 1;
gpx.option.jsn = tharg->option_jsn;