comparison Discovery/Src/buehlmann.c @ 224:ceecabfddb57 div-fixes-3

Bugfix, deco: fix 2 (small) problems with calculated ceiling This fixes 1 trivial, and 1 not really trivial bug in the calculation of the ceiling. When simulating a bounce dive to 80m, things become clear (tried this on a CCR dive, fixed setpoint 1.2bar, about 15 minutes of bottom time). Closely watch the behavior of the ceiling data. At some point during the ascent, the ceiling begins to decrease in 10cm steps. Then suddenly (while still ascending), the ceiling increases again with 1m, does not change for some time, and then suddenly steps 1.1m less deep. While not very relevant to real deco diving, it is simply wrong. The reason for this is subtle. The algorithm used to find the ceiling is a sort of linear search, stepping down a meter, overshoot the depth, and search back in 10cm steps. It seems some numerical instability. Fixing this, was a bit more computational intensive search by stepping up down in equal steps of 10cm. But, I'm pretty sure that things can be speeded up here, as a ceiling does not change fast, so it should be not that difficult to limit the search space, or use a binary search algorithm instead. The trivial second problem fixed, is that the ceiling ends at the surface and not at 1m depth. This small issue became visible after changing the step down size above. Signed-off-by: Jan Mulder <jlmulder@xs4all.nl>
author Jan Mulder <jlmulder@xs4all.nl>
date Sun, 31 Mar 2019 19:35:51 +0200
parents b7689d9e888a
children ff0d23625cd5
comparison
equal deleted inserted replaced
223:5f535ef6a3db 224:ceecabfddb57
13 #include <math.h> 13 #include <math.h>
14 #include <stdbool.h> 14 #include <stdbool.h>
15 #include "buehlmann.h" 15 #include "buehlmann.h"
16 #include "decom.h" 16 #include "decom.h"
17 17
18
19 extern const float helium_time_constant[16];
20 extern const float nitrogen_time_constant[16];
21 18
22 extern const float buehlmann_N2_a[]; 19 extern const float buehlmann_N2_a[];
23 extern const float buehlmann_N2_b[]; 20 extern const float buehlmann_N2_b[];
24 extern const float buehlmann_He_a[]; 21 extern const float buehlmann_He_a[];
25 extern const float buehlmann_He_b[]; 22 extern const float buehlmann_He_b[];
725 void buehlmann_ceiling_calculator(SLifeData* pLifeData, SDiveSettings * pDiveSettings, SDecoinfo * pDecoInfo) 722 void buehlmann_ceiling_calculator(SLifeData* pLifeData, SDiveSettings * pDiveSettings, SDecoinfo * pDecoInfo)
726 { 723 {
727 float gf_low; 724 float gf_low;
728 float gf_high; 725 float gf_high;
729 float gf_delta; 726 float gf_delta;
730 int dv_gf_low_stop_meter; 727 float dv_gf_low_stop_meter;
731 _Bool test_result; 728 _Bool test_result;
732 float next_gf_value; 729 float next_gf_value;
733 float next_pressure_absolute; 730 float next_pressure_absolute;
734 int depth_in_meter; 731 float depth_in_meter;
735 732
736 gf_low = pDiveSettings->gf_low; 733 gf_low = pDiveSettings->gf_low;
737 gf_high = pDiveSettings->gf_high; 734 gf_high = pDiveSettings->gf_high;
738 735
739 //
740 dv_gf_low_stop_meter = (int)((pDiveSettings->internal__pressure_first_stop_ambient_bar_as_upper_limit_for_gf_low_otherwise_zero - pLifeData->pressure_surface_bar) * 10); 736 dv_gf_low_stop_meter = (int)((pDiveSettings->internal__pressure_first_stop_ambient_bar_as_upper_limit_for_gf_low_otherwise_zero - pLifeData->pressure_surface_bar) * 10);
741 // 737
742 if(dv_gf_low_stop_meter < 1) 738 if(dv_gf_low_stop_meter < 1)
743 { 739 {
744 next_gf_value = gf_high; // fix hw 161024 740 next_gf_value = gf_high; // fix hw 161024
745 gf_delta = 0; 741 gf_delta = 0;
746 } 742 }
747 else 743 else
748 { 744 {
749 next_gf_value = gf_high; 745 next_gf_value = gf_high;
750 gf_delta = gf_high - gf_low; 746 gf_delta = gf_high - gf_low;
751 gf_delta /= dv_gf_low_stop_meter; // gf_delta is delta for each meter now!! 747 gf_delta /= (dv_gf_low_stop_meter * 10); // gf_delta is delta for 10 cm !!
752 } 748 }
753 // 749
754 depth_in_meter = 0; 750 depth_in_meter = 0;
755 next_pressure_absolute = pLifeData->pressure_surface_bar; 751 next_pressure_absolute = pLifeData->pressure_surface_bar;
756 752
757 memcpy(gTissue_nitrogen_bar, pLifeData->tissue_nitrogen_bar, (4*16)); 753 memcpy(gTissue_nitrogen_bar, pLifeData->tissue_nitrogen_bar, (4*16));
758 memcpy(gTissue_helium_bar, pLifeData->tissue_helium_bar, (4*16)); 754 memcpy(gTissue_helium_bar, pLifeData->tissue_helium_bar, (4*16));
760 // 756 //
761 test_result = buehlmann_tissue_test_tolerance(next_pressure_absolute); 757 test_result = buehlmann_tissue_test_tolerance(next_pressure_absolute);
762 // 758 //
763 while(!test_result && depth_in_meter < 200) 759 while(!test_result && depth_in_meter < 200)
764 { 760 {
765 depth_in_meter += 1; 761 depth_in_meter += 0.1;
766 next_gf_value = fmaxf(gf_low, next_gf_value - gf_delta); 762 next_gf_value = fmaxf(gf_low, next_gf_value - gf_delta);
767 gGF_value = next_gf_value / 100.0f; 763 gGF_value = next_gf_value / 100.0f;
768 next_pressure_absolute += 0.1f; // 1 meter down 764 next_pressure_absolute += 0.01f; // 0.1 meter down
769 test_result = buehlmann_tissue_test_tolerance(next_pressure_absolute); 765 test_result = buehlmann_tissue_test_tolerance(next_pressure_absolute);
770 } 766 }
771 // 767
772 if(test_result) 768 if(test_result)
773 { 769 {
774 // old direct paste
775 pDecoInfo->output_ceiling_meter = depth_in_meter; 770 pDecoInfo->output_ceiling_meter = depth_in_meter;
776 // new sub-meter hw 160331 771
777 if(depth_in_meter >= 1) 772 if(depth_in_meter >= 0)
778 { 773 {
779 for(int i = 0; i < 10; i++) 774 for(int i = 0; i < 10; i++)
780 { 775 {
781 next_gf_value += gf_delta/10.0f; 776 next_gf_value += gf_delta/10.0f;
782 gGF_value = next_gf_value / 100.0f; 777 gGF_value = next_gf_value / 100.0f;