HG_REPOSITORIES/LPP/INSTRUMENTATION/USERS/LEROY/LFR_basic-parameters Commit - r18:0f2eb26d750b

1

// In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)

1

// In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)

2

// version 1.0: 31/07/2013

2

// version 1.0: 31/07/2013

3

// version 1.1: 02/04/2014

3

// version 1.1: 02/04/2014

4

// version 1.2: 30/04/2014

4

// version 1.2: 30/04/2014

5

// version 1.3: 02/05/2014

5

// version 1.3: 02/05/2014

6

// version 1.4: 16/05/2014

6

// version 1.4: 16/05/2014

7

// version 1.5: 20/05/2014

7

// version 1.5: 20/05/2014

8

// version 1.6: 19/12/2014

8

// version 1.6: 19/12/2014

9

// version 1.7: 15/01/2015 (modifs de Paul + correction erreurs qui se compensaient (LSB <=> MSB + indices [0,2] <=> [1,3])

9

// version 1.7: 15/01/2015 (modifs de Paul + correction erreurs qui se compensaient (LSB <=> MSB + indices [0,2] <=> [1,3])

10

11

12

#include <stdint.h>

12

#include <stdint.h>

13

#include "basic_parameters.h"

13

#include "basic_parameters.h"

14

#include "basic_parameters_params.h"

14

#include "basic_parameters_params.h"

15

16

void init_k_coefficients(float *k_coefficients,

16

void init_k_coefficients(float *k_coefficients,

17

unsigned char nb_binscompressed_matrix )

17

unsigned char nb_binscompressed_matrix )

18

{

18

{

19

uint16_t i; // 16 bits unsigned

19

uint16_t i; // 16 bits unsigned

20

for(i=0; i<nb_binscompressed_matrix; i++){

20

for(i=0; i<nb_binscompressed_matrix; i++){

21

k_coefficients[i*NB_K_COEFF_PER_BIN+K44_PE] = 1;

21

k_coefficients[i*NB_K_COEFF_PER_BIN+K44_PE] = 1;

22

k_coefficients[i*NB_K_COEFF_PER_BIN+K55_PE] = 1;

22

k_coefficients[i*NB_K_COEFF_PER_BIN+K55_PE] = 1;

23

k_coefficients[i*NB_K_COEFF_PER_BIN+K45_PE_RE] = 1;

23

k_coefficients[i*NB_K_COEFF_PER_BIN+K45_PE_RE] = 1;

24

k_coefficients[i*NB_K_COEFF_PER_BIN+K45_PE_IM] = 1;

24

k_coefficients[i*NB_K_COEFF_PER_BIN+K45_PE_IM] = 1;

25

k_coefficients[i*NB_K_COEFF_PER_BIN+K14_SX_RE] = 1;

25

k_coefficients[i*NB_K_COEFF_PER_BIN+K14_SX_RE] = 1;

26

k_coefficients[i*NB_K_COEFF_PER_BIN+K14_SX_IM] = 1;

26

k_coefficients[i*NB_K_COEFF_PER_BIN+K14_SX_IM] = 1;

27

k_coefficients[i*NB_K_COEFF_PER_BIN+K15_SX_RE] = 1;

27

k_coefficients[i*NB_K_COEFF_PER_BIN+K15_SX_RE] = 1;

28

k_coefficients[i*NB_K_COEFF_PER_BIN+K15_SX_IM] = 1;

28

k_coefficients[i*NB_K_COEFF_PER_BIN+K15_SX_IM] = 1;

29

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_SX_RE] = 1;

29

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_SX_RE] = 1;

30

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_SX_IM] = 1;

30

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_SX_IM] = 1;

31

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_SX_RE] = 1;

31

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_SX_RE] = 1;

32

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_SX_IM] = 1;

32

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_SX_IM] = 1;

33

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_SX_RE] = 1;

33

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_SX_RE] = 1;

34

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_SX_IM] = 1;

34

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_SX_IM] = 1;

35

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_SX_RE] = 1;

35

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_SX_RE] = 1;

36

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_SX_IM] = 1;

36

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_SX_IM] = 1;

37

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NY_RE] = 1;

37

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NY_RE] = 1;

38

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NY_IM] = 1;

38

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NY_IM] = 1;

39

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NY_RE] = 1;

39

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NY_RE] = 1;

40

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NY_IM] = 1;

40

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NY_IM] = 1;

41

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NY_RE] = 1;

41

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NY_RE] = 1;

42

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NY_IM] = 1;

42

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NY_IM] = 1;

43

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NY_RE] = 1;

43

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NY_RE] = 1;

44

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NY_IM] = 1;

44

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NY_IM] = 1;

45

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] = 1;

45

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] = 1;

46

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] = 1;

46

k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] = 1;

47

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] = 1;

47

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] = 1;

48

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] = 1;

48

k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] = 1;

49

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] = 1;

49

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] = 1;

50

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] = 1;

50

k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] = 1;

51

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NZ_RE] = 1;

51

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NZ_RE] = 1;

52

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NZ_IM] = 1;

52

k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NZ_IM] = 1;

53

}

53

}

54

}

54

}

55

56

void BP1_set( float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1 ){

56

void BP1_set( float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1 ){

57

float PSDB; // 32-bit floating point

57

float PSDB; // 32-bit floating point

58

float PSDE;

58

float PSDE;

59

float tmp;

59

float tmp;

60

float NVEC_V0;

60

float NVEC_V0;

61

float NVEC_V1;

61

float NVEC_V1;

62

float NVEC_V2;

62

float NVEC_V2;

63

float aux;

63

float aux;

64

float tr_SB_SB;

64

float tr_SB_SB;

65

float e_cross_b_re;

65

float e_cross_b_re;

66

float e_cross_b_im;

66

float e_cross_b_im;

67

float n_cross_e_scal_b_re;

67

float n_cross_e_scal_b_re;

68

float n_cross_e_scal_b_im;

68

float n_cross_e_scal_b_im;

69

float ny;

69

float ny;

70

float nz;

70

float nz;

71

float bx_bx_star;

71

float bx_bx_star;

72

float vphi;

72

float vphi;

73

float significand;

73

float significand;

74

int exponent; // 32-bit signed integer

74

int exponent; // 32-bit signed integer

75

float alpha_M;

75

float alpha_M;

76

77

uint8_t nbitexp; // 8-bit unsigned integer

77

uint8_t nbitexp; // 8-bit unsigned integer

78

uint8_t nbitsig;

78

uint8_t nbitsig;

79

uint8_t tmp_uint8;

79

uint8_t tmp_uint8;

80

uint8_t *pt_uint8; // pointer on unsigned 8-bit integer

80

uint8_t *pt_uint8; // pointer on unsigned 8-bit integer

81

int8_t expmin; // 8-bit signed integer

81

int8_t expmin; // 8-bit signed integer

82

int8_t expmax;

82

int8_t expmax;

83

uint16_t rangesig; // 16-bit unsigned integer

83

uint16_t rangesig; // 16-bit unsigned integer

84

uint16_t psd;

84

uint16_t psd;

85

uint16_t exp;

85

uint16_t exp;

86

uint16_t tmp_uint16;

86

uint16_t tmp_uint16;

87

uint16_t i;

87

uint16_t i;

88

89

alpha_M = 45 * (3.1415927/180);

89

alpha_M = 45 * (3.1415927/180);

90

91

#ifdef DEBUG_TCH

91

#ifdef DEBUG_TCH

92

printf("BP1 : \n");

92

printf("BP1 : \n");

93

printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);

93

printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);

94

#endif

94

#endif

95

96

// initialization for managing the exponents of the floating point data:

96

// initialization for managing the exponents of the floating point data:

97

nbitexp = 5; // number of bits for the exponent

97

nbitexp = 5; // number of bits for the exponent

98

expmax = 30; // maximum value of the exponent

98

expmax = 30; // maximum value of the exponent

99

expmin = expmax - (1 << nbitexp) + 1; // accordingly the minimum exponent value

99

expmin = expmax - (1 << nbitexp) + 1; // accordingly the minimum exponent value

100

// for floating point data to be recorded on 12-bit words:

100

// for floating point data to be recorded on 12-bit words:

101

nbitsig = 12 - nbitexp; // number of bits for the significand

101

nbitsig = 12 - nbitexp; // number of bits for the significand

102

rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1

102

rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1

103

104

#ifdef DEBUG_TCH

104

#ifdef DEBUG_TCH

105

printf("nbitexp : %d, expmax : %d, expmin : %d\n", nbitexp, expmax, expmin);

105

printf("nbitexp : %d, expmax : %d, expmin : %d\n", nbitexp, expmax, expmin);

106

printf("nbitsig : %d, rangesig : %d\n", nbitsig, rangesig);

106

printf("nbitsig : %d, rangesig : %d\n", nbitsig, rangesig);

107

#endif

107

#endif

108

109

for(i=0; i<nb_bins_compressed_spec_mat; i++){

109

for(i=0; i<nb_bins_compressed_spec_mat; i++){

110

//==============================================

110

//==============================================

111

// BP1 PSDB == PA_LFR_SC_BP1_PB_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)

111

// BP1 PSDB == PA_LFR_SC_BP1_PB_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)

112

PSDB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] // S11

112

PSDB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] // S11

113

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22

113

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22

114

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]; // S33

114

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]; // S33

115

116

significand = frexpf(PSDB/3, &exponent); // 0.5 <= significand < 1

116

significand = frexpf(PSDB/3, &exponent); // 0.5 <= significand < 1

117

// PSDB/3 = significand * 2^exponent

117

// PSDB/3 = significand * 2^exponent

118

// the division by 3 is to ensure that max value <= 2^30

118

// the division by 3 is to ensure that max value <= 2^30

119

120

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

120

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

121

exponent = expmin;

121

exponent = expmin;

122

significand = 0.5; // min value that can be recorded

122

significand = 0.5; // min value that can be recorded

123

}

123

}

124

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

124

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

125

exponent = expmax;

125

exponent = expmax;

126

significand = 1.0; // max value that can be recorded

126

significand = 1.0; // max value that can be recorded

127

}

127

}

128

if (significand == 0) { // in that case exponent == 0 too

128

if (significand == 0) { // in that case exponent == 0 too

129

exponent = expmin;

129

exponent = expmin;

130

significand = 0.5; // min value that can be recorded

130

significand = 0.5; // min value that can be recorded

131

}

131

}

132

133

psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

133

psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

134

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

134

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

135

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

135

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

136

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

136

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

137

tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

137

tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

138

// left place of the significand bits (nbitsig),

138

// left place of the significand bits (nbitsig),

139

// making the 16-bit word to be recorded

139

// making the 16-bit word to be recorded

140

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

140

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

141

#ifdef MSB_FIRST_TCH

141

#ifdef MSB_FIRST_TCH

142

lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[0]; // Record MSB of tmp_uint16

142

lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[0]; // Record MSB of tmp_uint16

143

lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[1]; // Record LSB of tmp_uint16

143

lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[1]; // Record LSB of tmp_uint16

144

#endif

144

#endif

145

#ifdef LSB_FIRST_TCH

145

#ifdef LSB_FIRST_TCH

146

lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[1]; // Record MSB of tmp_uint16

146

lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[1]; // Record MSB of tmp_uint16

147

lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[0]; // Record LSB of tmp_uint16

147

lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[0]; // Record LSB of tmp_uint16

148

#endif

148

#endif

149

#ifdef DEBUG_TCH

149

#ifdef DEBUG_TCH

150

printf("\nBin number: %d\n", i);

150

printf("\nBin number: %d\n", i);

151

printf("PSDB / 3 : %16.8e\n",PSDB/3);

151

printf("PSDB / 3 : %16.8e\n",PSDB/3);

152

printf("significand : %16.8e\n",significand);

152

printf("significand : %16.8e\n",significand);

153

printf("exponent : %d\n" ,exponent);

153

printf("exponent : %d\n" ,exponent);

154

printf("psd for PSDB significand : %d\n",psd);

154

printf("psd for PSDB significand : %d\n",psd);

155

printf("exp for PSDB exponent : %d\n",exp);

155

printf("exp for PSDB exponent : %d\n",exp);

156

printf("pt_uint8[1] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]);

156

printf("pt_uint8[1] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]);

157

printf("pt_uint8[0] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]);

157

printf("pt_uint8[0] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]);

158

printf("lfr_bp1[i*NB_BYTES_BP1+3] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+3], lfr_bp1[i*NB_BYTES_BP1+3]);

158

printf("lfr_bp1[i*NB_BYTES_BP1+3] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+3], lfr_bp1[i*NB_BYTES_BP1+3]);

159

printf("lfr_bp1[i*NB_BYTES_BP1+2] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+2], lfr_bp1[i*NB_BYTES_BP1+2]);

159

printf("lfr_bp1[i*NB_BYTES_BP1+2] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+2], lfr_bp1[i*NB_BYTES_BP1+2]);

160

#endif

160

#endif

161

//==============================================

161

//==============================================

162

// BP1 PSDE == PA_LFR_SC_BP1_PE_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)

162

// BP1 PSDE == PA_LFR_SC_BP1_PE_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)

163

PSDE = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K44_PE] // S44

163

PSDE = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K44_PE] // S44

164

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K55_PE] // S55

164

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K55_PE] // S55

165

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K45_PE_RE] // S45 Re

165

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K45_PE_RE] // S45 Re

166

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K45_PE_IM]; // S45 Im

166

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K45_PE_IM]; // S45 Im

167

168

significand = frexpf(PSDE/2, &exponent); // 0.5 <= significand < 1

168

significand = frexpf(PSDE/2, &exponent); // 0.5 <= significand < 1

169

// PSDE/2 = significand * 2^exponent

169

// PSDE/2 = significand * 2^exponent

170

// the division by 2 is to ensure that max value <= 2^30

170

// the division by 2 is to ensure that max value <= 2^30

171

// should be reconsidered by taking into account the k-coefficients ...

171

// should be reconsidered by taking into account the k-coefficients ...

172

173

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

173

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

174

exponent = expmin;

174

exponent = expmin;

175

significand = 0.5; // min value that can be recorded

175

significand = 0.5; // min value that can be recorded

176

}

176

}

177

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

177

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

178

exponent = expmax;

178

exponent = expmax;

179

significand = 1.0; // max value that can be recorded

179

significand = 1.0; // max value that can be recorded

180

}

180

}

181

if (significand == 0) {// in that case exponent == 0 too

181

if (significand == 0) {// in that case exponent == 0 too

182

exponent = expmin;

182

exponent = expmin;

183

significand = 0.5; // min value that can be recorded

183

significand = 0.5; // min value that can be recorded

184

}

184

}

185

186

psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

186

psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

187

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

187

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

188

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

188

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

189

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

189

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

190

tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

190

tmp_uint16 = psd | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

191

// left place of the significand bits (nbitsig),

191

// left place of the significand bits (nbitsig),

192

// making the 16-bit word to be recorded

192

// making the 16-bit word to be recorded

193

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

193

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

194

#ifdef MSB_FIRST_TCH

194

#ifdef MSB_FIRST_TCH

195

lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[0]; // Record MSB of tmp_uint16

195

lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[0]; // Record MSB of tmp_uint16

196

lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[1]; // Record LSB of tmp_uint16

196

lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[1]; // Record LSB of tmp_uint16

197

#endif

197

#endif

198

#ifdef LSB_FIRST_TCH

198

#ifdef LSB_FIRST_TCH

199

lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[1]; // Record MSB of tmp_uint16

199

lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[1]; // Record MSB of tmp_uint16

200

lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[0]; // Record LSB of tmp_uint16

200

lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[0]; // Record LSB of tmp_uint16

201

#endif

201

#endif

202

#ifdef DEBUG_TCH

202

#ifdef DEBUG_TCH

203

printf("Bin number: %d\n", i);

203

printf("Bin number: %d\n", i);

204

printf("PSDE/2 : %16.8e\n",PSDE/2);

204

printf("PSDE/2 : %16.8e\n",PSDE/2);

205

printf("significand : %16.8e\n",significand);

205

printf("significand : %16.8e\n",significand);

206

printf("exponent : %d\n" ,exponent);

206

printf("exponent : %d\n" ,exponent);

207

printf("psd for PSDE significand : %d\n",psd);

207

printf("psd for PSDE significand : %d\n",psd);

208

printf("exp for PSDE exponent : %d\n",exp);

208

printf("exp for PSDE exponent : %d\n",exp);

209

printf("pt_uint8[1] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]);

209

printf("pt_uint8[1] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]);

210

printf("pt_uint8[0] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]);

210

printf("pt_uint8[0] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]);

211

printf("lfr_bp1[i*NB_BYTES_BP1+1] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+1], lfr_bp1[i*NB_BYTES_BP1+1]);

211

printf("lfr_bp1[i*NB_BYTES_BP1+1] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+1], lfr_bp1[i*NB_BYTES_BP1+1]);

212

printf("lfr_bp1[i*NB_BYTES_BP1+0] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+0], lfr_bp1[i*NB_BYTES_BP1+0]);

212

printf("lfr_bp1[i*NB_BYTES_BP1+0] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+0], lfr_bp1[i*NB_BYTES_BP1+0]);

213

#endif

213

#endif

214

//==============================================================================

214

//==============================================================================

215

// BP1 normal wave vector == PA_LFR_SC_BP1_NVEC_V0_F0 == 8 bits

215

// BP1 normal wave vector == PA_LFR_SC_BP1_NVEC_V0_F0 == 8 bits

216

// == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits

216

// == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits

217

// == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit

217

// == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit

218

tmp = sqrt( compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] //Im S12

218

tmp = sqrt_ple( compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] //Im S12

219

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] //Im S13

219

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] //Im S13

220

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] //Im S23

220

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] //Im S23

221

);

221

);

222

NVEC_V0 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]/ tmp; // S23 Im => n1

222

NVEC_V0 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]/ tmp; // S23 Im => n1

223

NVEC_V1 = -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / tmp; // S13 Im => n2

223

NVEC_V1 = -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / tmp; // S13 Im => n2

224

NVEC_V2 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / tmp; // S12 Im => n3

224

NVEC_V2 = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / tmp; // S12 Im => n3

225

226

lfr_bp1[i*NB_BYTES_BP1+4] = (uint8_t) (NVEC_V0*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

226

lfr_bp1[i*NB_BYTES_BP1+4] = (uint8_t) (NVEC_V0*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

227

lfr_bp1[i*NB_BYTES_BP1+5] = (uint8_t) (NVEC_V1*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

227

lfr_bp1[i*NB_BYTES_BP1+5] = (uint8_t) (NVEC_V1*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

228

pt_uint8 = (uint8_t*) &NVEC_V2; // Affect an uint8_t pointer with the adress of NVEC_V2

228

pt_uint8 = (uint8_t*) &NVEC_V2; // Affect an uint8_t pointer with the adress of NVEC_V2

229

#ifdef LSB_FIRST_TCH

229

#ifdef LSB_FIRST_TCH

230

lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[3] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention)

230

lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[3] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention)

231

// Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]

231

// Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]

232

#endif

232

#endif

233

#ifdef MSB_FIRST_TCH

233

#ifdef MSB_FIRST_TCH

234

lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[0] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 1th octet:SPARC convention)

234

lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[0] & 0x80; // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 1th octet:SPARC convention)

235

// Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]

235

// Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]

236

#endif

236

#endif

237

#ifdef DEBUG_TCH

237

#ifdef DEBUG_TCH

238

printf("NVEC_V0 : %16.8e\n",NVEC_V0);

238

printf("NVEC_V0 : %16.8e\n",NVEC_V0);

239

printf("NVEC_V1 : %16.8e\n",NVEC_V1);

239

printf("NVEC_V1 : %16.8e\n",NVEC_V1);

240

printf("NVEC_V2 : %16.8e\n",NVEC_V2);

240

printf("NVEC_V2 : %16.8e\n",NVEC_V2);

241

printf("lfr_bp1[i*NB_BYTES_BP1+4] for NVEC_V0 : %u\n",lfr_bp1[i*NB_BYTES_BP1+4]);

241

printf("lfr_bp1[i*NB_BYTES_BP1+4] for NVEC_V0 : %u\n",lfr_bp1[i*NB_BYTES_BP1+4]);

242

printf("lfr_bp1[i*NB_BYTES_BP1+5] for NVEC_V1 : %u\n",lfr_bp1[i*NB_BYTES_BP1+5]);

242

printf("lfr_bp1[i*NB_BYTES_BP1+5] for NVEC_V1 : %u\n",lfr_bp1[i*NB_BYTES_BP1+5]);

243

printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);

243

printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);

244

#endif

244

#endif

245

//=======================================================

245

//=======================================================

246

// BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits

246

// BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits

247

aux = 2*tmp / PSDB; // Compute the ellipticity

247

aux = 2*tmp / PSDB; // Compute the ellipticity

248

249

tmp_uint8 = (uint8_t) (aux*15 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

249

tmp_uint8 = (uint8_t) (aux*15 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

250

// where just the first 4 bits are used (0, ..., 15)

250

// where just the first 4 bits are used (0, ..., 15)

251

lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | (tmp_uint8 << 3); // Put these 4 bits next to the right place

251

lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | (tmp_uint8 << 3); // Put these 4 bits next to the right place

252

// of the sign bit of NVEC_V2 (recorded

252

// of the sign bit of NVEC_V2 (recorded

253

// previously in lfr_bp1[i*NB_BYTES_BP1+6])

253

// previously in lfr_bp1[i*NB_BYTES_BP1+6])

254

#ifdef DEBUG_TCH

254

#ifdef DEBUG_TCH

255

printf("ellipticity : %16.8e\n",aux);

255

printf("ellipticity : %16.8e\n",aux);

256

printf("tmp_uint8 for ellipticity : %u\n",tmp_uint8);

256

printf("tmp_uint8 for ellipticity : %u\n",tmp_uint8);

257

printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);

257

printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);

258

#endif

258

#endif

259

//==============================================================

259

//==============================================================

260

// BP1 degree of polarization == PA_LFR_SC_BP1_DOP_F0 == 3 bits

260

// BP1 degree of polarization == PA_LFR_SC_BP1_DOP_F0 == 3 bits

261

tr_SB_SB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]

261

tr_SB_SB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]

262

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]

262

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]

263

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]

263

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]

264

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1]

264

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1]

265

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2]

265

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2]

266

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3]

266

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3]

267

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4]

267

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4]

268

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]

268

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]

269

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11];

269

+ 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11];

270

aux = PSDB*PSDB;

270

aux = PSDB*PSDB;

271

tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux ); // Compute the degree of polarisation

271

tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux ); // Compute the degree of polarisation

272

273

tmp_uint8 = (uint8_t) (tmp*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

273

tmp_uint8 = (uint8_t) (tmp*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

274

// where just the first 3 bits are used (0, ..., 7)

274

// where just the first 3 bits are used (0, ..., 7)

275

lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | tmp_uint8; // Record these 3 bits at the 3 first bit positions

275

lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | tmp_uint8; // Record these 3 bits at the 3 first bit positions

276

// (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]

276

// (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]

277

#ifdef DEBUG_TCH

277

#ifdef DEBUG_TCH

278

printf("DOP : %16.8e\n",tmp);

278

printf("DOP : %16.8e\n",tmp);

279

printf("tmp_uint8 for DOP : %u\n",tmp_uint8);

279

printf("tmp_uint8 for DOP : %u\n",tmp_uint8);

280

printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity + DOP : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);

280

printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity + DOP : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);

281

#endif

281

#endif

282

//=======================================================================================

282

//=======================================================================================

283

// BP1 X_SO-component of the Poynting flux == PA_LFR_SC_BP1_SX_F0 == 8 (+ 2) bits

283

// BP1 X_SO-component of the Poynting flux == PA_LFR_SC_BP1_SX_F0 == 8 (+ 2) bits

284

// = 5 bits (exponent) + 3 bits (significand)

284

// = 5 bits (exponent) + 3 bits (significand)

285

// + 1 sign bit + 1 argument bit (two sectors)

285

// + 1 sign bit + 1 argument bit (two sectors)

286

e_cross_b_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_RE] //S34 Re

286

e_cross_b_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_RE] //S34 Re

287

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_RE] //S35 Re

287

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_RE] //S35 Re

288

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_RE] //S14 Re

288

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_RE] //S14 Re

289

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_RE] //S15 Re

289

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_RE] //S15 Re

290

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_RE] //S24 Re

290

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_RE] //S24 Re

291

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_RE] //S25 Re

291

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_RE] //S25 Re

292

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_IM] //S34 Im

292

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_IM] //S34 Im

293

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_IM] //S35 Im

293

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_IM] //S35 Im

294

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_IM] //S14 Im

294

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_IM] //S14 Im

295

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_IM] //S15 Im

295

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_IM] //S15 Im

296

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_IM] //S24 Im

296

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_IM] //S24 Im

297

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_IM]; //S25 Im

297

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_IM]; //S25 Im

298

// Im(S_ji) = -Im(S_ij)

298

// Im(S_ji) = -Im(S_ij)

299

// k_ji = k_ij

299

// k_ji = k_ij

300

e_cross_b_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_IM] //S34 Re

300

e_cross_b_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_IM] //S34 Re

301

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_IM] //S35 Re

301

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_IM] //S35 Re

302

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_IM] //S14 Re

302

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_IM] //S14 Re

303

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_IM] //S15 Re

303

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_IM] //S15 Re

304

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_IM] //S24 Re

304

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_IM] //S24 Re

305

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_IM] //S25 Re

305

+ compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_IM] //S25 Re

306

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_RE] //S34 Im

306

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_RE] //S34 Im

307

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_RE] //S35 Im

307

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_RE] //S35 Im

308

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_RE] //S14 Im

308

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_RE] //S14 Im

309

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_RE] //S15 Im

309

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_RE] //S15 Im

310

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_RE] //S24 Im

310

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_RE] //S24 Im

311

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_RE]; //S25 Im

311

- compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_RE]; //S25 Im

312

#ifdef DEBUG_TCH

312

#ifdef DEBUG_TCH

313

printf("ReaSX / 2 : %16.8e\n",e_cross_b_re/2);

313

printf("ReaSX / 2 : %16.8e\n",e_cross_b_re/2);

314

#endif

314

#endif

315

pt_uint8 = (uint8_t*) &e_cross_b_re; // Affect an uint8_t pointer with the adress of e_cross_b_re

315

pt_uint8 = (uint8_t*) &e_cross_b_re; // Affect an uint8_t pointer with the adress of e_cross_b_re

316

#ifdef LSB_FIRST_TCH

316

#ifdef LSB_FIRST_TCH

317

lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)

317

lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)

318

// Record it at the 8th bit position (from the right to the left)

318

// Record it at the 8th bit position (from the right to the left)

319

// of lfr_bp1[i*NB_BYTES_BP1+0]

319

// of lfr_bp1[i*NB_BYTES_BP1+0]

320

pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX|

320

pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX|

321

#endif

321

#endif

322

#ifdef MSB_FIRST_TCH

322

#ifdef MSB_FIRST_TCH

323

lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 1th octet:SPARC convention)

323

lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 1th octet:SPARC convention)

324

// Record it at the 8th bit position (from the right to the left)

324

// Record it at the 8th bit position (from the right to the left)

325

// of lfr_bp1[i*NB_BYTES_BP1+0]

325

// of lfr_bp1[i*NB_BYTES_BP1+0]

326

pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX|

326

pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make e_cross_b_re be positive in any case: |ReaSX|

327

#endif

327

#endif

328

significand = frexpf(e_cross_b_re/2, &exponent); // 0.5 <= significand < 1

328

significand = frexpf(e_cross_b_re/2, &exponent); // 0.5 <= significand < 1

329

// ReaSX/2 = significand * 2^exponent

329

// ReaSX/2 = significand * 2^exponent

330

// The division by 2 is to ensure that max value <= 2^30 (rough estimate)

330

// The division by 2 is to ensure that max value <= 2^30 (rough estimate)

331

// Should be reconsidered by taking into account the k-coefficients ...

331

// Should be reconsidered by taking into account the k-coefficients ...

332

333

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

333

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

334

exponent = expmin;

334

exponent = expmin;

335

significand = 0.5; // min value that can be recorded

335

significand = 0.5; // min value that can be recorded

336

}

336

}

337

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

337

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

338

exponent = expmax;

338

exponent = expmax;

339

significand = 1.0; // max value that can be recorded

339

significand = 1.0; // max value that can be recorded

340

}

340

}

341

if (significand == 0) { // in that case exponent == 0 too

341

if (significand == 0) { // in that case exponent == 0 too

342

exponent = expmin;

342

exponent = expmin;

343

significand = 0.5; // min value that can be recorded

343

significand = 0.5; // min value that can be recorded

344

}

344

}

345

346

lfr_bp1[i*NB_BYTES_BP1+7] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

346

lfr_bp1[i*NB_BYTES_BP1+7] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

347

// where just the first 3 bits are used (0, ..., 7)

347

// where just the first 3 bits are used (0, ..., 7)

348

tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where

348

tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where

349

// just the first 5 bits are used (0, ..., 2^5-1)

349

// just the first 5 bits are used (0, ..., 2^5-1)

350

#ifdef DEBUG_TCH

350

#ifdef DEBUG_TCH

351

printf("|ReaSX| / 2 : %16.8e\n",e_cross_b_re/2);

351

printf("|ReaSX| / 2 : %16.8e\n",e_cross_b_re/2);

352

printf("significand : %16.8e\n",significand);

352

printf("significand : %16.8e\n",significand);

353

printf("exponent : %d\n" ,exponent);

353

printf("exponent : %d\n" ,exponent);

354

printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]);

354

printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]);

355

printf("tmp_uint8 for ReaSX exponent : %d\n",tmp_uint8);

355

printf("tmp_uint8 for ReaSX exponent : %d\n",tmp_uint8);

356

#endif

356

#endif

357

lfr_bp1[i*NB_BYTES_BP1+7] = lfr_bp1[i*NB_BYTES_BP1+7] | (tmp_uint8 << 3); // Shift these 5 bits to the left before logical addition

357

lfr_bp1[i*NB_BYTES_BP1+7] = lfr_bp1[i*NB_BYTES_BP1+7] | (tmp_uint8 << 3); // Shift these 5 bits to the left before logical addition

358

// with lfr_bp1[i*NB_BYTES_BP1+7]

358

// with lfr_bp1[i*NB_BYTES_BP1+7]

359

#ifdef DEBUG_TCH

359

#ifdef DEBUG_TCH

360

printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]);

360

printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]);

361

printf("lfr_bp1[i*NB_BYTES_BP1+0] for ReaSX sign + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+0]);

361

printf("lfr_bp1[i*NB_BYTES_BP1+0] for ReaSX sign + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+0]);

362

printf("ImaSX / 2 : %16.8e\n",e_cross_b_im/2);

362

printf("ImaSX / 2 : %16.8e\n",e_cross_b_im/2);

363

#endif

363

#endif

364

pt_uint8 = (uint8_t*) &e_cross_b_im; // Affect an uint8_t pointer with the adress of e_cross_b_im

364

pt_uint8 = (uint8_t*) &e_cross_b_im; // Affect an uint8_t pointer with the adress of e_cross_b_im

365

#ifdef LSB_FIRST_TCH

365

#ifdef LSB_FIRST_TCH

366

pt_uint8[3] = pt_uint8[3] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 4th octet:PC convention)

366

pt_uint8[3] = pt_uint8[3] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 4th octet:PC convention)

367

#endif

367

#endif

368

#ifdef MSB_FIRST_TCH

368

#ifdef MSB_FIRST_TCH

369

pt_uint8[0] = pt_uint8[0] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 1th octet:SPARC convention)

369

pt_uint8[0] = pt_uint8[0] & 0x7f; // Make e_cross_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 1th octet:SPARC convention)

370

#endif

370

#endif

371

tmp_uint8 = (e_cross_b_im > e_cross_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect

371

tmp_uint8 = (e_cross_b_im > e_cross_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect

372

// an unsigned 8-bit char with 01000000; otherwise with null.

372

// an unsigned 8-bit char with 01000000; otherwise with null.

373

lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right

373

lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right

374

// to the left) of lfr_bp1[i*NB_BYTES_BP1+0], by simple logical addition.

374

// to the left) of lfr_bp1[i*NB_BYTES_BP1+0], by simple logical addition.

375

#ifdef DEBUG_TCH

375

#ifdef DEBUG_TCH

376

printf("|ImaSX| / 2 : %16.8e\n",e_cross_b_im/2);

376

printf("|ImaSX| / 2 : %16.8e\n",e_cross_b_im/2);

377

printf("ArgSX sign : %u\n",tmp_uint8);

377

printf("ArgSX sign : %u\n",tmp_uint8);

378

printf("lfr_bp1[i*NB_BYTES_BP1+0] for ReaSX & ArgSX signs + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+0]);

378

printf("lfr_bp1[i*NB_BYTES_BP1+0] for ReaSX & ArgSX signs + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+0]);

379

#endif

379

#endif

380

//======================================================================

380

//======================================================================

381

// BP1 phase velocity estimator == PA_LFR_SC_BP1_VPHI_F0 == 8 (+ 2) bits

381

// BP1 phase velocity estimator == PA_LFR_SC_BP1_VPHI_F0 == 8 (+ 2) bits

382

// = 5 bits (exponent) + 3 bits (significand)

382

// = 5 bits (exponent) + 3 bits (significand)

383

// + 1 sign bit + 1 argument bit (two sectors)

383

// + 1 sign bit + 1 argument bit (two sectors)

384

ny = sin(alpha_M)*NVEC_V1 + cos(alpha_M)*NVEC_V2;

384

ny = sin(alpha_M)*NVEC_V1 + cos(alpha_M)*NVEC_V2;

385

nz = NVEC_V0;

385

nz = NVEC_V0;

386

bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22 Re

386

bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9] // S22 Re

387

+ sin(alpha_M)*sin(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] // S33 Re

387

+ sin(alpha_M)*sin(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16] // S33 Re

388

- 2*sin(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]; // S23 Re

388

- 2*sin(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]; // S23 Re

389

390

n_cross_e_scal_b_re = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_RE] //S24 Re

390

n_cross_e_scal_b_re = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_RE] //S24 Re

391

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_RE] //S25 Re

391

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_RE] //S25 Re

392

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_RE] //S34 Re

392

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_RE] //S34 Re

393

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_RE] //S35 Re

393

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_RE] //S35 Re

394

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_IM] //S24 Im

394

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_IM] //S24 Im

395

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_IM] //S25 Im

395

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_IM] //S25 Im

396

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_IM] //S34 Im

396

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_IM] //S34 Im

397

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_IM]) //S35 Im

397

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_IM]) //S35 Im

398

+ nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] //S24 Re

398

+ nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] //S24 Re

399

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] //S25 Re

399

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] //S25 Re

400

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] //S34 Re

400

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] //S34 Re

401

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_RE] //S35 Re

401

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_RE] //S35 Re

402

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] //S24 Im

402

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] //S24 Im

403

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] //S25 Im

403

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] //S25 Im

404

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] //S34 Im

404

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] //S34 Im

405

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_IM]);//S35 Im

405

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_IM]);//S35 Im

406

// Im(S_ji) = -Im(S_ij)

406

// Im(S_ji) = -Im(S_ij)

407

// k_ji = k_ij

407

// k_ji = k_ij

408

n_cross_e_scal_b_im = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_IM] //S24 Re

408

n_cross_e_scal_b_im = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_IM] //S24 Re

409

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_IM] //S25 Re

409

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_IM] //S25 Re

410

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_IM] //S34 Re

410

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_IM] //S34 Re

411

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_IM] //S35 Re

411

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_IM] //S35 Re

412

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_RE] //S24 Im

412

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_RE] //S24 Im

413

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_RE] //S25 Im

413

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_RE] //S25 Im

414

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_RE] //S34 Im

414

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_RE] //S34 Im

415

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_RE]) //S35 Im

415

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_RE]) //S35 Im

416

+ nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] //S24 Re

416

+ nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] //S24 Re

417

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] //S25 Re

417

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] //S25 Re

418

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] //S34 Re

418

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] //S34 Re

419

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_IM] //S35 Re

419

+compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_IM] //S35 Re

420

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] //S24 Im

420

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] //S24 Im

421

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] //S25 Im

421

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] //S25 Im

422

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] //S34 Im

422

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] //S34 Im

423

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_RE]);//S35 Im

423

-compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_RE]);//S35 Im

424

#ifdef DEBUG_TCH

424

#ifdef DEBUG_TCH

425

printf("n_cross_e_scal_b_re : %16.8e\n",n_cross_e_scal_b_re);

425

printf("n_cross_e_scal_b_re : %16.8e\n",n_cross_e_scal_b_re);

426

printf("n_cross_e_scal_b_im : %16.8e\n",n_cross_e_scal_b_im);

426

printf("n_cross_e_scal_b_im : %16.8e\n",n_cross_e_scal_b_im);

427

#endif

427

#endif

428

// vphi = n_cross_e_scal_b_re / bx_bx_star => sign(VPHI) = sign(n_cross_e_scal_b_re)

428

// vphi = n_cross_e_scal_b_re / bx_bx_star => sign(VPHI) = sign(n_cross_e_scal_b_re)

429

pt_uint8 = (uint8_t*) &n_cross_e_scal_b_re; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_re

429

pt_uint8 = (uint8_t*) &n_cross_e_scal_b_re; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_re

430

#ifdef LSB_FIRST_TCH

430

#ifdef LSB_FIRST_TCH

431

lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)

431

lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | (pt_uint8[3] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)

432

// Record it at the 8th bit position (from the right to the left)

432

// Record it at the 8th bit position (from the right to the left)

433

// of lfr_bp1[i*NB_BYTES_BP1+2]

433

// of lfr_bp1[i*NB_BYTES_BP1+2]

434

pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|

434

pt_uint8[3] = (pt_uint8[3] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|

435

#endif

435

#endif

436

#ifdef MSB_FIRST_TCH

436

#ifdef MSB_FIRST_TCH

437

lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 1th octet:SPARC convention)

437

lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | (pt_uint8[0] & 0x80); // Extract its sign bit (32-bit float, sign bit in the 1th octet:SPARC convention)

438

// Record it at the 8th bit position (from the right to the left)

438

// Record it at the 8th bit position (from the right to the left)

439

// of lfr_bp1[i*NB_BYTES_BP1+2]

439

// of lfr_bp1[i*NB_BYTES_BP1+2]

440

pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|

440

pt_uint8[0] = (pt_uint8[0] & 0x7f); // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|

441

#endif

441

#endif

442

vphi = n_cross_e_scal_b_re / bx_bx_star; // Compute |VPHI|

442

vphi = n_cross_e_scal_b_re / bx_bx_star; // Compute |VPHI|

443

444

significand = frexpf(vphi/2, &exponent); // 0.5 <= significand < 1

444

significand = frexpf(vphi/2, &exponent); // 0.5 <= significand < 1

445

// vphi/2 = significand * 2^exponent

445

// vphi/2 = significand * 2^exponent

446

// The division by 2 is to ensure that max value <= 2^30 (rough estimate)

446

// The division by 2 is to ensure that max value <= 2^30 (rough estimate)

447

// Should be reconsidered by taking into account the k-coefficients ...

447

// Should be reconsidered by taking into account the k-coefficients ...

448

449

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

449

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

450

exponent = expmin;

450

exponent = expmin;

451

significand = 0.5; // min value that can be recorded

451

significand = 0.5; // min value that can be recorded

452

}

452

}

453

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

453

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

454

exponent = expmax;

454

exponent = expmax;

455

significand = 1.0; // max value that can be recorded

455

significand = 1.0; // max value that can be recorded

456

}

456

}

457

if (significand == 0) {// in that case exponent == 0 too

457

if (significand == 0) {// in that case exponent == 0 too

458

exponent = expmin;

458

exponent = expmin;

459

significand = 0.5; // min value that can be recorded

459

significand = 0.5; // min value that can be recorded

460

}

460

}

461

#ifdef DEBUG_TCH

461

#ifdef DEBUG_TCH

462

printf("|VPHI| / 2 : %16.8e\n",vphi/2);

462

printf("|VPHI| / 2 : %16.8e\n",vphi/2);

463

printf("significand : %16.8e\n",significand);

463

printf("significand : %16.8e\n",significand);

464

printf("exponent : %d\n" ,exponent);

464

printf("exponent : %d\n" ,exponent);

465

#endif

465

#endif

466

lfr_bp1[i*NB_BYTES_BP1+8] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

466

lfr_bp1[i*NB_BYTES_BP1+8] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding

467

// where just the first 3 bits are used (0, ..., 7)

467

// where just the first 3 bits are used (0, ..., 7)

468

tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where

468

tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where

469

// just the first 5 bits are used (0, ..., 2^5-1)

469

// just the first 5 bits are used (0, ..., 2^5-1)

470

#ifdef DEBUG_TCH

470

#ifdef DEBUG_TCH

471

printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]);

471

printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]);

472

printf("tmp_uint8 for VPHI exponent : %d\n",tmp_uint8);

472

printf("tmp_uint8 for VPHI exponent : %d\n",tmp_uint8);

473

#endif

473

#endif

474

lfr_bp1[i*NB_BYTES_BP1+8] = lfr_bp1[i*NB_BYTES_BP1+8] | (tmp_uint8 << 3); // shift these 5 bits to the left before logical addition

474

lfr_bp1[i*NB_BYTES_BP1+8] = lfr_bp1[i*NB_BYTES_BP1+8] | (tmp_uint8 << 3); // shift these 5 bits to the left before logical addition

475

// with lfr_bp1[i*NB_BYTES_BP1+8]

475

// with lfr_bp1[i*NB_BYTES_BP1+8]

476

#ifdef DEBUG_TCH

476

#ifdef DEBUG_TCH

477

printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]);

477

printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]);

478

printf("lfr_bp1[i*NB_BYTES_BP1+2] for VPHI sign + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+2]);

478

printf("lfr_bp1[i*NB_BYTES_BP1+2] for VPHI sign + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+2]);

479

#endif

479

#endif

480

pt_uint8 = (uint8_t*) &n_cross_e_scal_b_im; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_im

480

pt_uint8 = (uint8_t*) &n_cross_e_scal_b_im; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_im

481

#ifdef LSB_FIRST_TCH

481

#ifdef LSB_FIRST_TCH

482

pt_uint8[3] = pt_uint8[3] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 4th octet:PC convention)

482

pt_uint8[3] = pt_uint8[3] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 4th octet:PC convention)

483

#endif

483

#endif

484

#ifdef MSB_FIRST_TCH

484

#ifdef MSB_FIRST_TCH

485

pt_uint8[0] = pt_uint8[0] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 1th octet:SPARC convention)

485

pt_uint8[0] = pt_uint8[0] & 0x7f; // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 1th octet:SPARC convention)

486

#endif

486

#endif

487

tmp_uint8 = (n_cross_e_scal_b_im > n_cross_e_scal_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect

487

tmp_uint8 = (n_cross_e_scal_b_im > n_cross_e_scal_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect

488

// an unsigned 8-bit char with 01000000; otherwise with null.

488

// an unsigned 8-bit char with 01000000; otherwise with null.

489

lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right

489

lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right

490

// to the left) of lfr_bp1[i*NB_BYTES_BP1+3], by simple logical addition.

490

// to the left) of lfr_bp1[i*NB_BYTES_BP1+3], by simple logical addition.

491

#ifdef DEBUG_TCH

491

#ifdef DEBUG_TCH

492

printf("|n_cross_e_scal_b_im| : %16.8e\n",n_cross_e_scal_b_im);

492

printf("|n_cross_e_scal_b_im| : %16.8e\n",n_cross_e_scal_b_im);

493

printf("|n_cross_e_scal_b_im|/bx_bx_star/2: %16.8e\n",n_cross_e_scal_b_im/bx_bx_star/2);

493

printf("|n_cross_e_scal_b_im|/bx_bx_star/2: %16.8e\n",n_cross_e_scal_b_im/bx_bx_star/2);

494

printf("ArgNEBX sign : %u\n",tmp_uint8);

494

printf("ArgNEBX sign : %u\n",tmp_uint8);

495

printf("lfr_bp1[i*NB_BYTES_BP1+2] for VPHI & ArgNEBX signs + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+2]);

495

printf("lfr_bp1[i*NB_BYTES_BP1+2] for VPHI & ArgNEBX signs + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+2]);

496

#endif

496

#endif

497

}

497

}

498

}

498

}

499

500

float sqrt_ple( float x )

501

{

502

float y;

503

504

y = 0;

505

506

if ( x < 0)

507

{

508

printf("sqrt(<0)\n");

509

}

510

else

511

{

512

y = sqrt( x );

513

}

514

515

if ( y == 0 )

516

{

517

y = 1e-9;

518

}

519

520

return y;

521

}

522

500

void BP2_set( float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2 )

523

void BP2_set( float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2 )

501

{

524

{

502

float cross_re; // 32-bit floating point

525

float cross_re; // 32-bit floating point

503

float cross_im;

526

float cross_im;

504

float aux;

527

float aux;

505

float significand;

528

float significand;

506

int exponent; // 32-bit signed integer

529

int exponent; // 32-bit signed integer

507

uint8_t nbitexp; // 8-bit unsigned integer

530

uint8_t nbitexp; // 8-bit unsigned integer

508

uint8_t nbitsig;

531

uint8_t nbitsig;

509

uint8_t *pt_uint8; // pointer on unsigned 8-bit integer

532

uint8_t *pt_uint8; // pointer on unsigned 8-bit integer

510

int8_t expmin; // 8-bit signed integer

533

int8_t expmin; // 8-bit signed integer

511

int8_t expmax;

534

int8_t expmax;

512

uint16_t rangesig; // 16-bit unsigned integer

535

uint16_t rangesig; // 16-bit unsigned integer

513

uint16_t autocor;

536

uint16_t autocor;

514

uint16_t exp;

537

uint16_t exp;

515

uint16_t tmp_uint16;

538

uint16_t tmp_uint16;

516

uint16_t i;

539

uint16_t i;

517

540

518

#ifdef DEBUG_TCH

541

#ifdef DEBUG_TCH

519

printf("BP2 : \n");

542

printf("BP2 : \n");

520

printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);

543

printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);

521

#endif

544

#endif

522

545

523

// For floating point data to be recorded on 16-bit words :

546

// For floating point data to be recorded on 16-bit words :

524

nbitexp = 6; // number of bits for the exponent

547

nbitexp = 6; // number of bits for the exponent

525

nbitsig = 16 - nbitexp; // number of bits for the significand

548

nbitsig = 16 - nbitexp; // number of bits for the significand

526

rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1

549

rangesig = (1 << nbitsig)-1; // == 2^nbitsig - 1

527

expmax = 32;

550

expmax = 32;

528

expmin = expmax - (1 << nbitexp) + 1;

551

expmin = expmax - (1 << nbitexp) + 1;

529

552

530

#ifdef DEBUG_TCH

553

#ifdef DEBUG_TCH

531

printf("nbitexp : %d, nbitsig : %d, rangesig : %d\n", nbitexp, nbitsig, rangesig);

554

printf("nbitexp : %d, nbitsig : %d, rangesig : %d\n", nbitexp, nbitsig, rangesig);

532

printf("expmin : %d, expmax : %d\n", expmin, expmax);

555

printf("expmin : %d, expmax : %d\n", expmin, expmax);

533

#endif

556

#endif

534

557

535

for(i = 0; i<nb_bins_compressed_spec_mat; i++){

558

for(i = 0; i<nb_bins_compressed_spec_mat; i++){

536

//==============================================

559

//==============================================

537

// BP2 normalized cross correlations == PA_LFR_SC_BP2_CROSS_F0 == 10 * (8+8) bits

560

// BP2 normalized cross correlations == PA_LFR_SC_BP2_CROSS_F0 == 10 * (8+8) bits

538

// == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits

561

// == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits

539

// == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits

562

// == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits

540

// == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits

563

// == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits

541

// == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits

564

// == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits

542

// == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits

565

// == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits

543

// == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits

566

// == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits

544

// == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits

567

// == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits

545

// == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits

568

// == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits

546

// == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits

569

// == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits

547

// == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits

570

// == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits

548

// == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits

571

// == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits

549

// == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits

572

// == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits

550

// == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits

573

// == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits

551

// == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits

574

// == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits

552

// == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits

575

// == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits

553

// == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits

576

// == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits

554

// == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits

577

// == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits

555

// == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits

578

// == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits

556

// == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits

579

// == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits

557

// == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits

580

// == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits

558

// S12

581

// S12

559

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);

582

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);

560

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] / aux;

583

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] / aux;

561

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / aux;

584

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / aux;

562

lfr_bp2[i*NB_BYTES_BP2+10] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

585

lfr_bp2[i*NB_BYTES_BP2+10] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

563

lfr_bp2[i*NB_BYTES_BP2+20] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

586

lfr_bp2[i*NB_BYTES_BP2+20] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

564

#ifdef DEBUG_TCH

587

#ifdef DEBUG_TCH

565

printf("\nBin number: %d\n", i);

588

printf("\nBin number: %d\n", i);

566

printf("lfr_bp2[i*NB_BYTES_BP2+10] for cross12_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+10]);

589

printf("lfr_bp2[i*NB_BYTES_BP2+10] for cross12_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+10]);

567

printf("lfr_bp2[i*NB_BYTES_BP2+20] for cross12_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+20]);

590

printf("lfr_bp2[i*NB_BYTES_BP2+20] for cross12_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+20]);

568

#endif

591

#endif

569

// S13

592

// S13

570

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

593

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

571

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] / aux;

594

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] / aux;

572

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / aux;

595

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / aux;

573

lfr_bp2[i*NB_BYTES_BP2+11] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

596

lfr_bp2[i*NB_BYTES_BP2+11] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

574

lfr_bp2[i*NB_BYTES_BP2+21] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

597

lfr_bp2[i*NB_BYTES_BP2+21] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

575

#ifdef DEBUG_TCH

598

#ifdef DEBUG_TCH

576

printf("lfr_bp2[i*NB_BYTES_BP2+11] for cross13_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+11]);

599

printf("lfr_bp2[i*NB_BYTES_BP2+11] for cross13_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+11]);

577

printf("lfr_bp2[i*NB_BYTES_BP2+21] for cross13_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+21]);

600

printf("lfr_bp2[i*NB_BYTES_BP2+21] for cross13_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+21]);

578

#endif

601

#endif

579

// S14

602

// S14

580

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

603

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

581

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] / aux;

604

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] / aux;

582

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] / aux;

605

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] / aux;

583

lfr_bp2[i*NB_BYTES_BP2+12] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

606

lfr_bp2[i*NB_BYTES_BP2+12] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

584

lfr_bp2[i*NB_BYTES_BP2+22] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

607

lfr_bp2[i*NB_BYTES_BP2+22] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

585

#ifdef DEBUG_TCH

608

#ifdef DEBUG_TCH

586

printf("lfr_bp2[i*NB_BYTES_BP2+12] for cross14_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+12]);

609

printf("lfr_bp2[i*NB_BYTES_BP2+12] for cross14_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+12]);

587

printf("lfr_bp2[i*NB_BYTES_BP2+22] for cross14_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+22]);

610

printf("lfr_bp2[i*NB_BYTES_BP2+22] for cross14_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+22]);

588

#endif

611

#endif

589

// S15

612

// S15

590

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

613

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

591

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] / aux;

614

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] / aux;

592

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] / aux;

615

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] / aux;

593

lfr_bp2[i*NB_BYTES_BP2+13] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

616

lfr_bp2[i*NB_BYTES_BP2+13] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

594

lfr_bp2[i*NB_BYTES_BP2+23] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

617

lfr_bp2[i*NB_BYTES_BP2+23] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

595

#ifdef DEBUG_TCH

618

#ifdef DEBUG_TCH

596

printf("lfr_bp2[i*NB_BYTES_BP2+13] for cross15_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+13]);

619

printf("lfr_bp2[i*NB_BYTES_BP2+13] for cross15_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+13]);

597

printf("lfr_bp2[i*NB_BYTES_BP2+23] for cross15_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+23]);

620

printf("lfr_bp2[i*NB_BYTES_BP2+23] for cross15_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+23]);

598

#endif

621

#endif

599

// S23

622

// S23

600

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

623

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

601

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10] / aux;

624

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10] / aux;

602

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] / aux;

625

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] / aux;

603

lfr_bp2[i*NB_BYTES_BP2+14] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

626

lfr_bp2[i*NB_BYTES_BP2+14] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

604

lfr_bp2[i*NB_BYTES_BP2+24] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

627

lfr_bp2[i*NB_BYTES_BP2+24] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

605

#ifdef DEBUG_TCH

628

#ifdef DEBUG_TCH

606

printf("lfr_bp2[i*NB_BYTES_BP2+14] for cross23_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+14]);

629

printf("lfr_bp2[i*NB_BYTES_BP2+14] for cross23_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+14]);

607

printf("lfr_bp2[i*NB_BYTES_BP2+24] for cross23_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+24]);

630

printf("lfr_bp2[i*NB_BYTES_BP2+24] for cross23_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+24]);

608

#endif

631

#endif

609

// S24

632

// S24

610

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

633

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

611

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12] / aux;

634

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12] / aux;

612

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13] / aux;

635

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13] / aux;

613

lfr_bp2[i*NB_BYTES_BP2+15] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

636

lfr_bp2[i*NB_BYTES_BP2+15] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

614

lfr_bp2[i*NB_BYTES_BP2+25] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

637

lfr_bp2[i*NB_BYTES_BP2+25] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

615

#ifdef DEBUG_TCH

638

#ifdef DEBUG_TCH

616

printf("lfr_bp2[i*NB_BYTES_BP2+15] for cross24_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+15]);

639

printf("lfr_bp2[i*NB_BYTES_BP2+15] for cross24_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+15]);

617

printf("lfr_bp2[i*NB_BYTES_BP2+25] for cross24_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+25]);

640

printf("lfr_bp2[i*NB_BYTES_BP2+25] for cross24_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+25]);

618

#endif

641

#endif

619

// S25

642

// S25

620

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

643

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

621

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14] / aux;

644

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14] / aux;

622

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15] / aux;

645

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15] / aux;

623

lfr_bp2[i*NB_BYTES_BP2+16] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

646

lfr_bp2[i*NB_BYTES_BP2+16] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

624

lfr_bp2[i*NB_BYTES_BP2+26] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

647

lfr_bp2[i*NB_BYTES_BP2+26] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

625

#ifdef DEBUG_TCH

648

#ifdef DEBUG_TCH

626

printf("lfr_bp2[i*NB_BYTES_BP2+16] for cross25_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+16]);

649

printf("lfr_bp2[i*NB_BYTES_BP2+16] for cross25_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+16]);

627

printf("lfr_bp2[i*NB_BYTES_BP2+26] for cross25_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+26]);

650

printf("lfr_bp2[i*NB_BYTES_BP2+26] for cross25_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+26]);

628

#endif

651

#endif

629

// S34

652

// S34

630

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

653

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

631

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17] / aux;

654

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17] / aux;

632

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18] / aux;

655

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18] / aux;

633

lfr_bp2[i*NB_BYTES_BP2+17] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

656

lfr_bp2[i*NB_BYTES_BP2+17] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

634

lfr_bp2[i*NB_BYTES_BP2+27] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

657

lfr_bp2[i*NB_BYTES_BP2+27] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

635

#ifdef DEBUG_TCH

658

#ifdef DEBUG_TCH

636

printf("lfr_bp2[i*NB_BYTES_BP2+17] for cross34_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+17]);

659

printf("lfr_bp2[i*NB_BYTES_BP2+17] for cross34_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+17]);

637

printf("lfr_bp2[i*NB_BYTES_BP2+27] for cross34_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+27]);

660

printf("lfr_bp2[i*NB_BYTES_BP2+27] for cross34_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+27]);

638

#endif

661

#endif

639

// S35

662

// S35

640

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

663

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

641

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19] / aux;

664

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19] / aux;

642

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20] / aux;

665

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20] / aux;

643

lfr_bp2[i*NB_BYTES_BP2+18] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

666

lfr_bp2[i*NB_BYTES_BP2+18] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

644

lfr_bp2[i*NB_BYTES_BP2+28] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

667

lfr_bp2[i*NB_BYTES_BP2+28] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

645

#ifdef DEBUG_TCH

668

#ifdef DEBUG_TCH

646

printf("lfr_bp2[i*NB_BYTES_BP2+18] for cross35_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+18]);

669

printf("lfr_bp2[i*NB_BYTES_BP2+18] for cross35_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+18]);

647

printf("lfr_bp2[i*NB_BYTES_BP2+28] for cross35_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+28]);

670

printf("lfr_bp2[i*NB_BYTES_BP2+28] for cross35_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+28]);

648

#endif

671

#endif

649

// S45

672

// S45

650

aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

673

aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

651

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] / aux;

674

cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] / aux;

652

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] / aux;

675

cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] / aux;

653

lfr_bp2[i*NB_BYTES_BP2+19] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

676

lfr_bp2[i*NB_BYTES_BP2+19] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

654

lfr_bp2[i*NB_BYTES_BP2+29] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

677

lfr_bp2[i*NB_BYTES_BP2+29] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding

655

#ifdef DEBUG_TCH

678

#ifdef DEBUG_TCH

656

printf("lfr_bp2[i*NB_BYTES_BP2+19] for cross45_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+19]);

679

printf("lfr_bp2[i*NB_BYTES_BP2+19] for cross45_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+19]);

657

printf("lfr_bp2[i*NB_BYTES_BP2+29] for cross45_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+29]);

680

printf("lfr_bp2[i*NB_BYTES_BP2+29] for cross45_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+29]);

658

#endif

681

#endif

659

//==============================================

682

//==============================================

660

// BP2 auto correlations == PA_LFR_SC_BP2_AUTO_F0 == 5*16 bits = 5*[6 bits (exponent) + 10 bits (significand)]

683

// BP2 auto correlations == PA_LFR_SC_BP2_AUTO_F0 == 5*16 bits = 5*[6 bits (exponent) + 10 bits (significand)]

661

// == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits

684

// == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits

662

// == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits

685

// == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits

663

// == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits

686

// == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits

664

// == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits

687

// == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits

665

// == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits

688

// == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits

666

// S11

689

// S11

667

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX], &exponent); // 0.5 <= significand < 1

690

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX], &exponent); // 0.5 <= significand < 1

668

// S11 = significand * 2^exponent

691

// S11 = significand * 2^exponent

669

#ifdef DEBUG_TCH

692

#ifdef DEBUG_TCH

670

printf("S11 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]);

693

printf("S11 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]);

671

printf("significand : %16.8e\n",significand);

694

printf("significand : %16.8e\n",significand);

672

printf("exponent : %d\n" ,exponent);

695

printf("exponent : %d\n" ,exponent);

673

#endif

696

#endif

674

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

697

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

675

exponent = expmin;

698

exponent = expmin;

676

significand = 0.5; // min value that can be recorded

699

significand = 0.5; // min value that can be recorded

677

}

700

}

678

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

701

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

679

exponent = expmax;

702

exponent = expmax;

680

significand = 1.0; // max value that can be recorded

703

significand = 1.0; // max value that can be recorded

681

}

704

}

682

if (significand == 0) { // in that case exponent == 0 too

705

if (significand == 0) { // in that case exponent == 0 too

683

exponent = expmin;

706

exponent = expmin;

684

significand = 0.5; // min value that can be recorded

707

significand = 0.5; // min value that can be recorded

685

}

708

}

686

709

687

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

710

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

688

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

711

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

689

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

712

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

690

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

713

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

691

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

714

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

692

// left place of the significand bits (nbitsig),

715

// left place of the significand bits (nbitsig),

693

// making the 16-bit word to be recorded

716

// making the 16-bit word to be recorded

694

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

717

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

695

#ifdef MSB_FIRST_TCH

718

#ifdef MSB_FIRST_TCH

696

lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[0]; // Record MSB of tmp_uint16

719

lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[0]; // Record MSB of tmp_uint16

697

lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[1]; // Record LSB of tmp_uint16

720

lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[1]; // Record LSB of tmp_uint16

698

#endif

721

#endif

699

#ifdef LSB_FIRST_TCH

722

#ifdef LSB_FIRST_TCH

700

lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[1]; // Record MSB of tmp_uint16

723

lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[1]; // Record MSB of tmp_uint16

701

lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[0]; // Record LSB of tmp_uint16

724

lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[0]; // Record LSB of tmp_uint16

702

#endif

725

#endif

703

#ifdef DEBUG_TCH

726

#ifdef DEBUG_TCH

704

printf("autocor for S11 significand : %u\n",autocor);

727

printf("autocor for S11 significand : %u\n",autocor);

705

printf("exp for S11 exponent : %u\n",exp);

728

printf("exp for S11 exponent : %u\n",exp);

706

printf("pt_uint8[1] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

729

printf("pt_uint8[1] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

707

printf("pt_uint8[0] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

730

printf("pt_uint8[0] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

708

printf("lfr_bp2[i*NB_BYTES_BP2+1] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+1], lfr_bp2[i*NB_BYTES_BP2+1]);

731

printf("lfr_bp2[i*NB_BYTES_BP2+1] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+1], lfr_bp2[i*NB_BYTES_BP2+1]);

709

printf("lfr_bp2[i*NB_BYTES_BP2+0] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+0], lfr_bp2[i*NB_BYTES_BP2+0]);

732

printf("lfr_bp2[i*NB_BYTES_BP2+0] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+0], lfr_bp2[i*NB_BYTES_BP2+0]);

710

#endif

733

#endif

711

// S22

734

// S22

712

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9], &exponent); // 0.5 <= significand < 1

735

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9], &exponent); // 0.5 <= significand < 1

713

// S22 = significand * 2^exponent

736

// S22 = significand * 2^exponent

714

#ifdef DEBUG_TCH

737

#ifdef DEBUG_TCH

715

printf("S22 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);

738

printf("S22 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);

716

printf("significand : %16.8e\n",significand);

739

printf("significand : %16.8e\n",significand);

717

printf("exponent : %d\n" ,exponent);

740

printf("exponent : %d\n" ,exponent);

718

#endif

741

#endif

719

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

742

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

720

exponent = expmin;

743

exponent = expmin;

721

significand = 0.5; // min value that can be recorded

744

significand = 0.5; // min value that can be recorded

722

}

745

}

723

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

746

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

724

exponent = expmax;

747

exponent = expmax;

725

significand = 1.0; // max value that can be recorded

748

significand = 1.0; // max value that can be recorded

726

}

749

}

727

if (significand == 0) { // in that case exponent == 0 too

750

if (significand == 0) { // in that case exponent == 0 too

728

exponent = expmin;

751

exponent = expmin;

729

significand = 0.5; // min value that can be recorded

752

significand = 0.5; // min value that can be recorded

730

}

753

}

731

754

732

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

755

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

733

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

756

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

734

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

757

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

735

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

758

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

736

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

759

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

737

// left place of the significand bits (nbitsig),

760

// left place of the significand bits (nbitsig),

738

// making the 16-bit word to be recorded

761

// making the 16-bit word to be recorded

739

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

762

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

740

#ifdef MSB_FIRST_TCH

763

#ifdef MSB_FIRST_TCH

741

lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[0]; // Record MSB of tmp_uint16

764

lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[0]; // Record MSB of tmp_uint16

742

lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[1]; // Record LSB of tmp_uint16

765

lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[1]; // Record LSB of tmp_uint16

743

#endif

766

#endif

744

#ifdef LSB_FIRST_TCH

767

#ifdef LSB_FIRST_TCH

745

lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[1]; // Record MSB of tmp_uint16

768

lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[1]; // Record MSB of tmp_uint16

746

lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[0]; // Record LSB of tmp_uint16

769

lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[0]; // Record LSB of tmp_uint16

747

#endif

770

#endif

748

#ifdef DEBUG_TCH

771

#ifdef DEBUG_TCH

749

printf("autocor for S22 significand : %u\n",autocor);

772

printf("autocor for S22 significand : %u\n",autocor);

750

printf("exp for S11 exponent : %u\n",exp);

773

printf("exp for S11 exponent : %u\n",exp);

751

printf("pt_uint8[1] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

774

printf("pt_uint8[1] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

752

printf("pt_uint8[0] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

775

printf("pt_uint8[0] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

753

printf("lfr_bp2[i*NB_BYTES_BP2+3] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+3], lfr_bp2[i*NB_BYTES_BP2+3]);

776

printf("lfr_bp2[i*NB_BYTES_BP2+3] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+3], lfr_bp2[i*NB_BYTES_BP2+3]);

754

printf("lfr_bp2[i*NB_BYTES_BP2+2] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+2], lfr_bp2[i*NB_BYTES_BP2+2]);

777

printf("lfr_bp2[i*NB_BYTES_BP2+2] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+2], lfr_bp2[i*NB_BYTES_BP2+2]);

755

#endif

778

#endif

756

// S33

779

// S33

757

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16], &exponent); // 0.5 <= significand < 1

780

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16], &exponent); // 0.5 <= significand < 1

758

// S33 = significand * 2^exponent

781

// S33 = significand * 2^exponent

759

#ifdef DEBUG_TCH

782

#ifdef DEBUG_TCH

760

printf("S33 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

783

printf("S33 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

761

printf("significand : %16.8e\n",significand);

784

printf("significand : %16.8e\n",significand);

762

printf("exponent : %d\n" ,exponent);

785

printf("exponent : %d\n" ,exponent);

763

#endif

786

#endif

764

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

787

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

765

exponent = expmin;

788

exponent = expmin;

766

significand = 0.5; // min value that can be recorded

789

significand = 0.5; // min value that can be recorded

767

}

790

}

768

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

791

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

769

exponent = expmax;

792

exponent = expmax;

770

significand = 1.0; // max value that can be recorded

793

significand = 1.0; // max value that can be recorded

771

}

794

}

772

if (significand == 0) { // in that case exponent == 0 too

795

if (significand == 0) { // in that case exponent == 0 too

773

exponent = expmin;

796

exponent = expmin;

774

significand = 0.5; // min value that can be recorded

797

significand = 0.5; // min value that can be recorded

775

}

798

}

776

799

777

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

800

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

778

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

801

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

779

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

802

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

780

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

803

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

781

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

804

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

782

// left place of the significand bits (nbitsig),

805

// left place of the significand bits (nbitsig),

783

// making the 16-bit word to be recorded

806

// making the 16-bit word to be recorded

784

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

807

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

785

#ifdef MSB_FIRST_TCH

808

#ifdef MSB_FIRST_TCH

786

lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[0]; // Record MSB of tmp_uint16

809

lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[0]; // Record MSB of tmp_uint16

787

lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[1]; // Record LSB of tmp_uint16

810

lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[1]; // Record LSB of tmp_uint16

788

#endif

811

#endif

789

#ifdef LSB_FIRST_TCH

812

#ifdef LSB_FIRST_TCH

790

lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[1]; // Record MSB of tmp_uint16

813

lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[1]; // Record MSB of tmp_uint16

791

lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[0]; // Record LSB of tmp_uint16

814

lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[0]; // Record LSB of tmp_uint16

792

#endif

815

#endif

793

#ifdef DEBUG_TCH

816

#ifdef DEBUG_TCH

794

printf("autocor for S33 significand : %u\n",autocor);

817

printf("autocor for S33 significand : %u\n",autocor);

795

printf("exp for S33 exponent : %u\n",exp);

818

printf("exp for S33 exponent : %u\n",exp);

796

printf("pt_uint8[1] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

819

printf("pt_uint8[1] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

797

printf("pt_uint8[0] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

820

printf("pt_uint8[0] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

798

printf("lfr_bp2[i*NB_BYTES_BP2+5] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+5], lfr_bp2[i*NB_BYTES_BP2+5]);

821

printf("lfr_bp2[i*NB_BYTES_BP2+5] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+5], lfr_bp2[i*NB_BYTES_BP2+5]);

799

printf("lfr_bp2[i*NB_BYTES_BP2+4] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+4], lfr_bp2[i*NB_BYTES_BP2+4]);

822

printf("lfr_bp2[i*NB_BYTES_BP2+4] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+4], lfr_bp2[i*NB_BYTES_BP2+4]);

800

#endif

823

#endif

801

// S44

824

// S44

802

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21], &exponent); // 0.5 <= significand < 1

825

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21], &exponent); // 0.5 <= significand < 1

803

// S44 = significand * 2^exponent

826

// S44 = significand * 2^exponent

804

#ifdef DEBUG_TCH

827

#ifdef DEBUG_TCH

805

printf("S44 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

828

printf("S44 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);

806

printf("significand : %16.8e\n",significand);

829

printf("significand : %16.8e\n",significand);

807

printf("exponent : %d\n" ,exponent);

830

printf("exponent : %d\n" ,exponent);

808

#endif

831

#endif

809

832

810

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

833

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

811

exponent = expmin;

834

exponent = expmin;

812

significand = 0.5; // min value that can be recorded

835

significand = 0.5; // min value that can be recorded

813

}

836

}

814

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

837

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

815

exponent = expmax;

838

exponent = expmax;

816

significand = 1.0; // max value that can be recorded

839

significand = 1.0; // max value that can be recorded

817

}

840

}

818

if (significand == 0) { // in that case exponent == 0 too

841

if (significand == 0) { // in that case exponent == 0 too

819

exponent = expmin;

842

exponent = expmin;

820

significand = 0.5; // min value that can be recorded

843

significand = 0.5; // min value that can be recorded

821

}

844

}

822

845

823

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

846

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

824

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

847

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

825

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

848

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

826

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

849

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

827

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

850

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

828

// left place of the significand bits (nbitsig),

851

// left place of the significand bits (nbitsig),

829

// making the 16-bit word to be recorded

852

// making the 16-bit word to be recorded

830

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

853

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

831

#ifdef MSB_FIRST_TCH

854

#ifdef MSB_FIRST_TCH

832

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record MSB of tmp_uint16

855

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record MSB of tmp_uint16

833

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record LSB of tmp_uint16

856

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record LSB of tmp_uint16

834

#endif

857

#endif

835

#ifdef LSB_FIRST_TCH

858

#ifdef LSB_FIRST_TCH

836

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record MSB of tmp_uint16

859

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record MSB of tmp_uint16

837

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record LSB of tmp_uint16

860

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record LSB of tmp_uint16

838

#endif

861

#endif

839

#ifdef DEBUG_TCH

862

#ifdef DEBUG_TCH

840

printf("autocor for S44 significand : %u\n",autocor);

863

printf("autocor for S44 significand : %u\n",autocor);

841

printf("exp for S44 exponent : %u\n",exp);

864

printf("exp for S44 exponent : %u\n",exp);

842

printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

865

printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

843

printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

866

printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

844

printf("lfr_bp2[i*NB_BYTES_BP2+7] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+7], lfr_bp2[i*NB_BYTES_BP2+7]);

867

printf("lfr_bp2[i*NB_BYTES_BP2+7] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+7], lfr_bp2[i*NB_BYTES_BP2+7]);

845

printf("lfr_bp2[i*NB_BYTES_BP2+6] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+6], lfr_bp2[i*NB_BYTES_BP2+6]);

868

printf("lfr_bp2[i*NB_BYTES_BP2+6] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+6], lfr_bp2[i*NB_BYTES_BP2+6]);

846

#endif

869

#endif

847

// S55

870

// S55

848

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent); // 0.5 <= significand < 1

871

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent); // 0.5 <= significand < 1

849

// S55 = significand * 2^exponent

872

// S55 = significand * 2^exponent

850

#ifdef DEBUG_TCH

873

#ifdef DEBUG_TCH

851

printf("S55 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

874

printf("S55 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

852

printf("significand : %16.8e\n",significand);

875

printf("significand : %16.8e\n",significand);

853

printf("exponent : %d\n" ,exponent);

876

printf("exponent : %d\n" ,exponent);

854

#endif

877

#endif

855

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

878

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

856

exponent = expmin;

879

exponent = expmin;

857

significand = 0.5; // min value that can be recorded

880

significand = 0.5; // min value that can be recorded

858

}

881

}

859

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

882

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

860

exponent = expmax;

883

exponent = expmax;

861

significand = 1.0; // max value that can be recorded

884

significand = 1.0; // max value that can be recorded

862

}

885

}

863

if (significand == 0) { // in that case exponent == 0 too

886

if (significand == 0) { // in that case exponent == 0 too

864

exponent = expmin;

887

exponent = expmin;

865

significand = 0.5; // min value that can be recorded

888

significand = 0.5; // min value that can be recorded

866

}

889

}

867

890

868

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

891

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

869

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

892

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

870

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

893

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

871

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

894

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

872

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

895

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

873

// left place of the significand bits (nbitsig),

896

// left place of the significand bits (nbitsig),

874

// making the 16-bit word to be recorded

897

// making the 16-bit word to be recorded

875

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

898

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

876

#ifdef MSB_FIRST_TCH

899

#ifdef MSB_FIRST_TCH

877

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record MSB of tmp_uint16

900

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record MSB of tmp_uint16

878

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record LSB of tmp_uint16

901

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record LSB of tmp_uint16

879

//printf("MSB:\n");

902

//printf("MSB:\n");

880

#endif

903

#endif

881

#ifdef LSB_FIRST_TCH

904

#ifdef LSB_FIRST_TCH

882

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record MSB of tmp_uint16

905

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record MSB of tmp_uint16

883

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record LSB of tmp_uint16

906

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record LSB of tmp_uint16

884

//printf("LSB:\n");

907

//printf("LSB:\n");

885

#endif

908

#endif

886

#ifdef DEBUG_TCH

909

#ifdef DEBUG_TCH

887

printf("autocor for S55 significand : %u\n",autocor);

910

printf("autocor for S55 significand : %u\n",autocor);

888

printf("exp for S55 exponent : %u\n",exp);

911

printf("exp for S55 exponent : %u\n",exp);

889

printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

912

printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

890

printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

913

printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

891

printf("lfr_bp2[i*NB_BYTES_BP2+9] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+9], lfr_bp2[i*NB_BYTES_BP2+9]);

914

printf("lfr_bp2[i*NB_BYTES_BP2+9] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+9], lfr_bp2[i*NB_BYTES_BP2+9]);

892

printf("lfr_bp2[i*NB_BYTES_BP2+8] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+8], lfr_bp2[i*NB_BYTES_BP2+8]);

915

printf("lfr_bp2[i*NB_BYTES_BP2+8] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+8], lfr_bp2[i*NB_BYTES_BP2+8]);

893

#endif

916

#endif

894

}

917

}

895

}

918

}

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages

             // In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)
             // version 1.0: 31/07/2013
             // version 1.1: 02/04/2014
             // version 1.2: 30/04/2014
             // version 1.3: 02/05/2014
             // version 1.4: 16/05/2014
             // version 1.5: 20/05/2014
             // version 1.6: 19/12/2014
             // version 1.7: 15/01/2015 (modifs de Paul + correction erreurs qui se compensaient (LSB <=> MSB + indices [0,2] <=> [1,3])
             #include <stdint.h>
             #include "basic_parameters.h"
             #include "basic_parameters_params.h"
             void init_k_coefficients(float *k_coefficients,
                                      unsigned char nb_binscompressed_matrix )
             {
                 uint16_t i; // 16 bits unsigned
                 for(i=0; i<nb_binscompressed_matrix; i++){
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K44_PE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K55_PE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K45_PE_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K45_PE_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K14_SX_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K14_SX_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K15_SX_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K15_SX_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K24_SX_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K24_SX_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K25_SX_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K25_SX_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K34_SX_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K34_SX_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K35_SX_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K35_SX_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NY_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NY_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NY_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NY_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NY_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NY_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NY_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NY_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NZ_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K24_NZ_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NZ_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K25_NZ_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NZ_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K34_NZ_IM] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NZ_RE] = 1;
                     k_coefficients[i*NB_K_COEFF_PER_BIN+K35_NZ_IM] = 1;
                 }
             }
             void BP1_set( float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1 ){
                 float PSDB;                         // 32-bit floating point
                 float PSDE;
                 float tmp;
                 float NVEC_V0;
                 float NVEC_V1;
                 float NVEC_V2;
                 float aux;
                 float tr_SB_SB;
                 float e_cross_b_re;
                 float e_cross_b_im;
                 float n_cross_e_scal_b_re;
                 float n_cross_e_scal_b_im;
                 float ny;
                 float nz;
                 float bx_bx_star;
                 float vphi;
                 float significand;
                 int exponent;                       // 32-bit signed integer
                 float alpha_M;
                 uint8_t nbitexp;                    // 8-bit unsigned integer
                 uint8_t nbitsig;
                 uint8_t tmp_uint8;
                 uint8_t *pt_uint8;                  // pointer on unsigned 8-bit integer
                 int8_t expmin;                      // 8-bit signed integer
                 int8_t expmax;
                 uint16_t rangesig;                  // 16-bit unsigned integer
                 uint16_t psd;
                 uint16_t exp;
                 uint16_t tmp_uint16;
                 uint16_t i;
                 alpha_M = 45 * (3.1415927/180);
             #ifdef DEBUG_TCH
                 printf("BP1 : \n");
                 printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);
             #endif
                 // initialization for managing the exponents of the floating point data:
                 nbitexp = 5;                           // number of bits for the exponent
                 expmax = 30;                           // maximum value of the exponent
                 expmin = expmax - (1 << nbitexp) + 1;  // accordingly the minimum exponent value
                 // for floating point data to be recorded on 12-bit words:
                 nbitsig = 12 - nbitexp;       // number of bits for the significand
                 rangesig = (1 << nbitsig)-1;  // == 2^nbitsig - 1
             #ifdef DEBUG_TCH
                 printf("nbitexp : %d, expmax : %d, expmin : %d\n", nbitexp, expmax, expmin);
                 printf("nbitsig : %d, rangesig : %d\n", nbitsig, rangesig);
             #endif
                 for(i=0; i<nb_bins_compressed_spec_mat; i++){
                     //==============================================
                     // BP1 PSDB == PA_LFR_SC_BP1_PB_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)
                     PSDB =  compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]          // S11
                           + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]        // S22
                           + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16];      // S33
                     significand = frexpf(PSDB/3, &exponent);  // 0.5 <= significand < 1
                                                               // PSDB/3 = significand * 2^exponent
                     // the division by 3 is to ensure that max value <= 2^30
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                          // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = psd | (exp << nbitsig);     // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[1]; // Record LSB of tmp_uint16
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+2] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp1[i*NB_BYTES_BP1+3] = pt_uint8[0]; // Record LSB of tmp_uint16
             #endif
             #ifdef DEBUG_TCH
                     printf("\nBin number: %d\n", i);
                     printf("PSDB / 3    : %16.8e\n",PSDB/3);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
                     printf("psd for PSDB significand : %d\n",psd);
                     printf("exp for PSDB exponent : %d\n",exp);
                     printf("pt_uint8[1] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for PSDB exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+3] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+3], lfr_bp1[i*NB_BYTES_BP1+3]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+2] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+2], lfr_bp1[i*NB_BYTES_BP1+2]);
             #endif
                     //==============================================
                     // BP1 PSDE == PA_LFR_SC_BP1_PE_F0 == 12 bits = 5 bits (exponent) + 7 bits (significand)
                     PSDE =  compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K44_PE]        // S44
                           + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K55_PE]        // S55
                           + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K45_PE_RE]     // S45 Re
                           - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] * k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K45_PE_IM];    // S45 Im
                     significand = frexpf(PSDE/2, &exponent); // 0.5 <= significand < 1
                                                              // PSDE/2 = significand * 2^exponent
                     // the division by 2 is to ensure that max value <= 2^30
                     // should be reconsidered by taking into account the k-coefficients ...
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {// in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;   // min value that can be recorded
                     }
                     psd = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                          // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = psd | (exp << nbitsig);     // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[1]; // Record LSB of tmp_uint16
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+0] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp1[i*NB_BYTES_BP1+1] = pt_uint8[0]; // Record LSB of tmp_uint16
             #endif
             #ifdef DEBUG_TCH
                     printf("Bin number: %d\n", i);
                     printf("PSDE/2      : %16.8e\n",PSDE/2);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
                     printf("psd for PSDE significand : %d\n",psd);
                     printf("exp for PSDE exponent : %d\n",exp);
                     printf("pt_uint8[1] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for PSDE exponent + significand: %.3d or %.2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+1] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+1], lfr_bp1[i*NB_BYTES_BP1+1]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+0] : %.3d or %.2x\n",lfr_bp1[i*NB_BYTES_BP1+0], lfr_bp1[i*NB_BYTES_BP1+0]);
             #endif
                     //==============================================================================
                     // BP1 normal wave vector == PA_LFR_SC_BP1_NVEC_V0_F0 == 8 bits
                                            // == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits
                                            // == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit
-                    tmp = sqrt( compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2]   //Im S12
+                    tmp = sqrt_ple( compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2]   //Im S12
                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4]   //Im S13
                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]  //Im S23
                                );
                     NVEC_V0 =  compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]/ tmp;  // S23 Im  => n1
                     NVEC_V1 = -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / tmp;  // S13 Im  => n2
                     NVEC_V2 =  compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / tmp;  // S12 Im  => n3
                     lfr_bp1[i*NB_BYTES_BP1+4] = (uint8_t) (NVEC_V0*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp1[i*NB_BYTES_BP1+5] = (uint8_t) (NVEC_V1*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     pt_uint8 = (uint8_t*) &NVEC_V2;                              // Affect an uint8_t pointer with the adress of NVEC_V2
             #ifdef LSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[3] & 0x80;  // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 4th octet:PC convention)
                                                                      // Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]
             #endif
             #ifdef MSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+6] = pt_uint8[0] & 0x80;  // Extract the sign bit of NVEC_V2 (32-bit float, sign bit in the 1th octet:SPARC convention)
                                                                      // Record it at the 8th bit position (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]
             #endif
             #ifdef DEBUG_TCH
                     printf("NVEC_V0  : %16.8e\n",NVEC_V0);
                     printf("NVEC_V1  : %16.8e\n",NVEC_V1);
                     printf("NVEC_V2  : %16.8e\n",NVEC_V2);
                     printf("lfr_bp1[i*NB_BYTES_BP1+4] for NVEC_V0 : %u\n",lfr_bp1[i*NB_BYTES_BP1+4]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+5] for NVEC_V1 : %u\n",lfr_bp1[i*NB_BYTES_BP1+5]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);
             #endif
                     //=======================================================
                     // BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits
                     aux = 2*tmp / PSDB;                   // Compute the ellipticity
                     tmp_uint8 = (uint8_t) (aux*15 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding
                                                           // where just the first 4 bits are used (0, ..., 15)
                     lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | (tmp_uint8 << 3); // Put these 4 bits next to the right place
                                                                                               // of the sign bit of NVEC_V2 (recorded
                                                                                               // previously in lfr_bp1[i*NB_BYTES_BP1+6])
             #ifdef DEBUG_TCH
                     printf("ellipticity  : %16.8e\n",aux);
                     printf("tmp_uint8 for ellipticity : %u\n",tmp_uint8);
                     printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);
             #endif
                     //==============================================================
                     // BP1 degree of polarization == PA_LFR_SC_BP1_DOP_F0 == 3 bits
                     tr_SB_SB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]     *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]
                              + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]   *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]
                              + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]  *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]
                              + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1]
                              + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2]
                              + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3]
                              + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] *compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4]
                              + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]
                              + 2 * compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11];
                     aux = PSDB*PSDB;
                     tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux );  // Compute the degree of polarisation
                     tmp_uint8 = (uint8_t) (tmp*7 + 0.5);       // Shift and cast into a 8-bit uint8_t with rounding
                                                                // where just the first 3 bits are used (0, ..., 7)
                     lfr_bp1[i*NB_BYTES_BP1+6] = lfr_bp1[i*NB_BYTES_BP1+6] | tmp_uint8; // Record these 3 bits at the 3 first bit positions
                                                                                        // (from the right to the left) of lfr_bp1[i*NB_BYTES_BP1+6]
             #ifdef DEBUG_TCH
                     printf("DOP  : %16.8e\n",tmp);
                     printf("tmp_uint8 for DOP : %u\n",tmp_uint8);
                     printf("lfr_bp1[i*NB_BYTES_BP1+6] for NVEC_V2 + ellipticity + DOP : %u\n",lfr_bp1[i*NB_BYTES_BP1+6]);
             #endif
                     //=======================================================================================
                     // BP1 X_SO-component of the Poynting flux == PA_LFR_SC_BP1_SX_F0 == 8 (+ 2) bits
                     //                                          = 5 bits (exponent) + 3 bits (significand)
                     //                                          + 1 sign bit + 1 argument bit (two sectors)
                     e_cross_b_re =  compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_RE]  //S34 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_RE]  //S35 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_RE]  //S14 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_RE]  //S15 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_RE]  //S24 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_RE]  //S25 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_IM]  //S34 Im
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_IM]  //S35 Im
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_IM]  //S14 Im
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_IM]  //S15 Im
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_IM]  //S24 Im
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_IM]; //S25 Im
                     // Im(S_ji) = -Im(S_ij)
                     // k_ji = k_ij
                     e_cross_b_im =  compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_IM]  //S34 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_IM]  //S35 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_IM]  //S14 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_IM]  //S15 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_IM]  //S24 Re
                                   + compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_IM]  //S25 Re
                                   - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_SX_RE]  //S34 Im
                                   - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_SX_RE]  //S35 Im
                                   - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K14_SX_RE]  //S14 Im
                                   - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] *k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K15_SX_RE]  //S15 Im
                                   - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_SX_RE]  //S24 Im
                                   - compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_SX_RE]; //S25 Im
             #ifdef DEBUG_TCH
                     printf("ReaSX / 2   : %16.8e\n",e_cross_b_re/2);
             #endif
                     pt_uint8 = (uint8_t*) &e_cross_b_re;      // Affect an uint8_t pointer with the adress of e_cross_b_re
             #ifdef LSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | (pt_uint8[3] & 0x80);  // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)
                                                                                       // Record it at the 8th bit position (from the right to the left)
                                                                                       // of lfr_bp1[i*NB_BYTES_BP1+0]
                     pt_uint8[3] = (pt_uint8[3] & 0x7f);       // Make e_cross_b_re be positive in any case: |ReaSX|
             #endif
             #ifdef MSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] | (pt_uint8[0] & 0x80);  // Extract its sign bit (32-bit float, sign bit in the 1th octet:SPARC convention)
                                                                                       // Record it at the 8th bit position (from the right to the left)
                                                                                       // of lfr_bp1[i*NB_BYTES_BP1+0]
                     pt_uint8[0] = (pt_uint8[0] & 0x7f);       // Make e_cross_b_re be positive in any case: |ReaSX|
             #endif
                     significand = frexpf(e_cross_b_re/2, &exponent); // 0.5 <= significand < 1
                                                                      // ReaSX/2 = significand * 2^exponent
                     // The division by 2 is to ensure that max value <= 2^30 (rough estimate)
                     // Should be reconsidered by taking into account the k-coefficients ...
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     lfr_bp1[i*NB_BYTES_BP1+7] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding
                                                                                        // where just the first 3 bits are used (0, ..., 7)
                     tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where
                                                              // just the first 5 bits are used (0, ..., 2^5-1)
             #ifdef DEBUG_TCH
                     printf("|ReaSX| / 2 : %16.8e\n",e_cross_b_re/2);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
                     printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]);
                     printf("tmp_uint8        for ReaSX exponent    : %d\n",tmp_uint8);
             #endif
                     lfr_bp1[i*NB_BYTES_BP1+7] = lfr_bp1[i*NB_BYTES_BP1+7] | (tmp_uint8 << 3); // Shift these 5 bits to the left before logical addition
                                                                                               // with lfr_bp1[i*NB_BYTES_BP1+7]
             #ifdef DEBUG_TCH
                     printf("lfr_bp1[i*NB_BYTES_BP1+7] for ReaSX exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+7]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+0] for ReaSX sign + PSDE 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+0]);
                     printf("ImaSX / 2   : %16.8e\n",e_cross_b_im/2);
             #endif
                     pt_uint8 = (uint8_t*) &e_cross_b_im; // Affect an uint8_t pointer with the adress of e_cross_b_im
             #ifdef LSB_FIRST_TCH
                     pt_uint8[3] = pt_uint8[3] & 0x7f;    // Make e_cross_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 4th octet:PC convention)
             #endif
             #ifdef MSB_FIRST_TCH
                     pt_uint8[0] = pt_uint8[0] & 0x7f;    // Make e_cross_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 1th octet:SPARC convention)
             #endif
                     tmp_uint8 = (e_cross_b_im > e_cross_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect
                                                                              // an unsigned 8-bit char with 01000000; otherwise with null.
                     lfr_bp1[i*NB_BYTES_BP1+0] = lfr_bp1[i*NB_BYTES_BP1+0] |  tmp_uint8; // Record it as a sign bit at the 7th bit position (from the right
                                                                                         // to the left) of lfr_bp1[i*NB_BYTES_BP1+0], by simple logical addition.
             #ifdef DEBUG_TCH
                     printf("|ImaSX| / 2 : %16.8e\n",e_cross_b_im/2);
                     printf("ArgSX sign  : %u\n",tmp_uint8);
                     printf("lfr_bp1[i*NB_BYTES_BP1+0] for ReaSX & ArgSX signs + PSDE 'exponent'  : %u\n",lfr_bp1[i*NB_BYTES_BP1+0]);
             #endif
                     //======================================================================
                     // BP1 phase velocity estimator == PA_LFR_SC_BP1_VPHI_F0 == 8 (+ 2) bits
                     //                                          = 5 bits (exponent) + 3 bits (significand)
                     //                                          + 1 sign bit + 1 argument bit (two sectors)
                     ny = sin(alpha_M)*NVEC_V1 + cos(alpha_M)*NVEC_V2;
                     nz = NVEC_V0;
                     bx_bx_star = cos(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]   // S22 Re
                                + sin(alpha_M)*sin(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]  // S33 Re
                              - 2*sin(alpha_M)*cos(alpha_M)*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10]; // S23 Re
                     n_cross_e_scal_b_re = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_RE]  //S24 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_RE]  //S25 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_RE]  //S34 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_RE]  //S35 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_IM]  //S24 Im
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_IM]  //S25 Im
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_IM]  //S34 Im
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_IM]) //S35 Im
                                         + nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_RE]  //S24 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_RE]  //S25 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_RE]  //S34 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_RE]  //S35 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_IM]  //S24 Im
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_IM]  //S25 Im
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_IM]  //S34 Im
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_IM]);//S35 Im
                     // Im(S_ji) = -Im(S_ij)
                     // k_ji = k_ij
                     n_cross_e_scal_b_im = ny * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_IM]  //S24 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_IM]  //S25 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_IM]  //S34 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_IM]  //S35 Re
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NY_RE]  //S24 Im
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NY_RE]  //S25 Im
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NY_RE]  //S34 Im
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NY_RE]) //S35 Im
                                         + nz * (compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_IM]  //S24 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_IM]  //S25 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_IM]  //S34 Re
                                                +compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_IM]  //S35 Re
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K24_NZ_RE]  //S24 Im
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K25_NZ_RE]  //S25 Im
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K34_NZ_RE]  //S34 Im
                                                -compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20]*k_coeff_intercalib[i*NB_K_COEFF_PER_BIN+K35_NZ_RE]);//S35 Im
             #ifdef DEBUG_TCH
                     printf("n_cross_e_scal_b_re   : %16.8e\n",n_cross_e_scal_b_re);
                     printf("n_cross_e_scal_b_im   : %16.8e\n",n_cross_e_scal_b_im);
             #endif
                     // vphi = n_cross_e_scal_b_re / bx_bx_star => sign(VPHI) = sign(n_cross_e_scal_b_re)
                     pt_uint8 = (uint8_t*) &n_cross_e_scal_b_re; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_re
             #ifdef LSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | (pt_uint8[3] & 0x80);  // Extract its sign bit (32-bit float, sign bit in the 4th octet:PC convention)
                                                                                     // Record it at the 8th bit position (from the right to the left)
                                                                                     // of lfr_bp1[i*NB_BYTES_BP1+2]
                     pt_uint8[3] = (pt_uint8[3] & 0x7f);     // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|
             #endif
             #ifdef MSB_FIRST_TCH
                     lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] | (pt_uint8[0] & 0x80);  // Extract its sign bit (32-bit float, sign bit in the 1th octet:SPARC convention)
                                                                                     // Record it at the 8th bit position (from the right to the left)
                                                                                     // of lfr_bp1[i*NB_BYTES_BP1+2]
                     pt_uint8[0] = (pt_uint8[0] & 0x7f);     // Make n_cross_e_scal_b_re be positive in any case: |n_cross_e_scal_b_re|
             #endif
                     vphi = n_cross_e_scal_b_re / bx_bx_star;  // Compute |VPHI|
                     significand = frexpf(vphi/2, &exponent);  // 0.5 <= significand < 1
                                                               // vphi/2 = significand * 2^exponent
                     // The division by 2 is to ensure that max value <= 2^30 (rough estimate)
                     // Should be reconsidered by taking into account the k-coefficients ...
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {// in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;   // min value that can be recorded
                     }
             #ifdef DEBUG_TCH
                     printf("|VPHI| / 2  : %16.8e\n",vphi/2);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
             #endif
                     lfr_bp1[i*NB_BYTES_BP1+8] = (uint8_t) ((significand*2-1)*7 + 0.5); // Shift and cast into a 8-bit uint8_t with rounding
                                                                                        // where just the first 3 bits are used (0, ..., 7)
                     tmp_uint8 = (uint8_t) (exponent-expmin); // Shift and cast into a 8-bit uint8_t where
                                                              // just the first 5 bits are used (0, ..., 2^5-1)
             #ifdef DEBUG_TCH
                     printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]);
                     printf("tmp_uint8        for VPHI exponent    : %d\n",tmp_uint8);
             #endif
                     lfr_bp1[i*NB_BYTES_BP1+8] = lfr_bp1[i*NB_BYTES_BP1+8] | (tmp_uint8 << 3); // shift these 5 bits to the left before logical addition
                                                                                // with lfr_bp1[i*NB_BYTES_BP1+8]
             #ifdef DEBUG_TCH
                     printf("lfr_bp1[i*NB_BYTES_BP1+8] for VPHI exponent + significand : %u\n",lfr_bp1[i*NB_BYTES_BP1+8]);
                     printf("lfr_bp1[i*NB_BYTES_BP1+2] for VPHI sign + PSDB 'exponent' : %u\n",lfr_bp1[i*NB_BYTES_BP1+2]);
             #endif
                     pt_uint8 = (uint8_t*) &n_cross_e_scal_b_im; // Affect an uint8_t pointer with the adress of n_cross_e_scal_b_im
             #ifdef LSB_FIRST_TCH
                     pt_uint8[3] = pt_uint8[3] & 0x7f;           // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 4th octet:PC convention)
             #endif
             #ifdef MSB_FIRST_TCH
                     pt_uint8[0] = pt_uint8[0] & 0x7f;           // Make n_cross_e_scal_b_im be positive in any case: |ImaSX| (32-bit float, sign bit in the 1th octet:SPARC convention)
             #endif
                     tmp_uint8 = (n_cross_e_scal_b_im > n_cross_e_scal_b_re) ? 0x40 : 0x00; // Determine the sector argument of SX. If |Im| > |Re| affect
                                                                                            // an unsigned 8-bit char with 01000000; otherwise with null.
                     lfr_bp1[i*NB_BYTES_BP1+2] = lfr_bp1[i*NB_BYTES_BP1+2] |  tmp_uint8;    // Record it as a sign bit at the 7th bit position (from the right
                                                                                            // to the left) of lfr_bp1[i*NB_BYTES_BP1+3], by simple logical addition.
             #ifdef DEBUG_TCH
                     printf("|n_cross_e_scal_b_im|             : %16.8e\n",n_cross_e_scal_b_im);
                     printf("|n_cross_e_scal_b_im|/bx_bx_star/2: %16.8e\n",n_cross_e_scal_b_im/bx_bx_star/2);
                     printf("ArgNEBX sign          : %u\n",tmp_uint8);
                     printf("lfr_bp1[i*NB_BYTES_BP1+2] for VPHI & ArgNEBX signs + PSDB 'exponent'  : %u\n",lfr_bp1[i*NB_BYTES_BP1+2]);
             #endif
                 }
             }
+            float sqrt_ple( float x )
+            {
+                float y;
+                y = 0;
+                if ( x < 0)
+                {
+                    printf("sqrt(<0)\n");
+                }
+                else
+                {
+                    y = sqrt( x );
+                }
+                if ( y  == 0 )
+                {
+                    y = 1e-9;
+                }
+                return y;
+            }
             void BP2_set( float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2 )
             {
                 float cross_re;                     // 32-bit floating point
                 float cross_im;
                 float aux;
                 float significand;
                 int exponent;                       // 32-bit signed integer
                 uint8_t nbitexp;                    // 8-bit unsigned integer
                 uint8_t nbitsig;
                 uint8_t *pt_uint8;                  // pointer on unsigned 8-bit integer
                 int8_t expmin;                      // 8-bit signed integer
                 int8_t expmax;
                 uint16_t rangesig;                  // 16-bit unsigned integer
                 uint16_t autocor;
                 uint16_t exp;
                 uint16_t tmp_uint16;
                 uint16_t i;
             #ifdef DEBUG_TCH
                 printf("BP2 : \n");
                 printf("Number of bins: %d\n", nb_bins_compressed_spec_mat);
             #endif
                 // For floating point data to be recorded on 16-bit words :
                 nbitexp = 6;                  // number of bits for the exponent
                 nbitsig = 16 - nbitexp;       // number of bits for the significand
                 rangesig = (1 << nbitsig)-1;  // == 2^nbitsig - 1
                 expmax = 32;
                 expmin = expmax - (1 << nbitexp) + 1;
             #ifdef DEBUG_TCH
                 printf("nbitexp : %d, nbitsig : %d, rangesig : %d\n", nbitexp, nbitsig, rangesig);
                 printf("expmin : %d, expmax : %d\n", expmin, expmax);
             #endif
                 for(i = 0; i<nb_bins_compressed_spec_mat; i++){
                     //==============================================
                     // BP2 normalized cross correlations == PA_LFR_SC_BP2_CROSS_F0 == 10 * (8+8) bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits
                                                       // == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits
                     // S12
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+1] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+2] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+10] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+20] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("\nBin number: %d\n", i);
                     printf("lfr_bp2[i*NB_BYTES_BP2+10] for cross12_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+10]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+20] for cross12_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+20]);
             #endif
                     // S13
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+3] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+4] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+11] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+21] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+11] for cross13_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+11]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+21] for cross13_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+21]);
             #endif
                     // S14
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+5] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+6] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+12] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+22] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+12] for cross14_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+12]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+22] for cross14_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+22]);
             #endif
                     // S15
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+7] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+8] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+13] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+23] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+13] for cross15_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+13]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+23] for cross15_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+23]);
             #endif
                     // S23
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+10] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+11] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+14] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+24] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+14] for cross23_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+14]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+24] for cross23_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+24]);
             #endif
                     // S24
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+12] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+13] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+15] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+25] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+15] for cross24_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+15]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+25] for cross24_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+25]);
             #endif
                     // S25
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+14] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+15] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+16] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+26] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+16] for cross25_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+16]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+26] for cross25_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+26]);
             #endif
                     // S34
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+17] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+18] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+17] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+27] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+17] for cross34_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+17]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+27] for cross34_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+27]);
             #endif
                     // S35
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+19] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+20] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+18] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+28] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+18] for cross35_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+18]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+28] for cross35_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+28]);
             #endif
                     // S45
-                    aux = sqrt(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
+                    aux = sqrt_ple(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
                     cross_re = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+22] / aux;
                     cross_im = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+23] / aux;
                     lfr_bp2[i*NB_BYTES_BP2+19] = (uint8_t) (cross_re*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
                     lfr_bp2[i*NB_BYTES_BP2+29] = (uint8_t) (cross_im*127.5 + 128); // Shift and cast into a 8-bit uint8_t (0, ..., 255) with rounding
             #ifdef DEBUG_TCH
                     printf("lfr_bp2[i*NB_BYTES_BP2+19] for cross45_re (%16.8e) : %.3u\n",cross_re, lfr_bp2[i*NB_BYTES_BP2+19]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+29] for cross45_im (%16.8e) : %.3u\n",cross_im, lfr_bp2[i*NB_BYTES_BP2+29]);
             #endif
                     //==============================================
                     // BP2  auto correlations == PA_LFR_SC_BP2_AUTO_F0 == 5*16 bits = 5*[6 bits (exponent) + 10 bits (significand)]
                                            // == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits
                                            // == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits
                                            // == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits
                                            // == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits
                                            // == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits
                     // S11
                     significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX], &exponent);  // 0.5 <= significand < 1
                                                                                                             // S11 = significand * 2^exponent
             #ifdef DEBUG_TCH
                     printf("S11         : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX]);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
             #endif
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                              // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[1]; // Record LSB of tmp_uint16
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+0] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+1] = pt_uint8[0]; // Record LSB of tmp_uint16
             #endif
             #ifdef DEBUG_TCH
                     printf("autocor for S11 significand : %u\n",autocor);
                     printf("exp for S11 exponent : %u\n",exp);
                     printf("pt_uint8[1] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for S11 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+1] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+1], lfr_bp2[i*NB_BYTES_BP2+1]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+0] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+0], lfr_bp2[i*NB_BYTES_BP2+0]);
             #endif
                     // S22
                     significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9], &exponent);  // 0.5 <= significand < 1
                                                                                                               // S22 = significand * 2^exponent
             #ifdef DEBUG_TCH
                     printf("S22         : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
             #endif
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                              // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[1]; // Record LSB of tmp_uint16
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+2] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+3] = pt_uint8[0]; // Record LSB of tmp_uint16
             #endif
             #ifdef DEBUG_TCH
                     printf("autocor for S22 significand : %u\n",autocor);
                     printf("exp for S11 exponent : %u\n",exp);
                     printf("pt_uint8[1] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for S22 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+3] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+3], lfr_bp2[i*NB_BYTES_BP2+3]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+2] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+2], lfr_bp2[i*NB_BYTES_BP2+2]);
             #endif
                     // S33
                     significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16], &exponent);  // 0.5 <= significand < 1
                                                                                                                // S33 = significand * 2^exponent
             #ifdef DEBUG_TCH
                     printf("S33         : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
             #endif
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                              // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[1]; // Record LSB of tmp_uint16
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+4] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+5] = pt_uint8[0]; // Record LSB of tmp_uint16
             #endif
             #ifdef DEBUG_TCH
                     printf("autocor for S33 significand : %u\n",autocor);
                     printf("exp for S33 exponent : %u\n",exp);
                     printf("pt_uint8[1] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for S33 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+5] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+5], lfr_bp2[i*NB_BYTES_BP2+5]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+4] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+4], lfr_bp2[i*NB_BYTES_BP2+4]);
             #endif
                     // S44
                     significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21], &exponent);  // 0.5 <= significand < 1
                                                                                                                // S44 = significand * 2^exponent
             #ifdef DEBUG_TCH
                     printf("S44         : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+21]);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
             #endif
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                              // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record LSB of tmp_uint16
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record LSB of tmp_uint16
             #endif
             #ifdef DEBUG_TCH
                     printf("autocor for S44 significand : %u\n",autocor);
                     printf("exp for S44 exponent : %u\n",exp);
                     printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+7] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+7], lfr_bp2[i*NB_BYTES_BP2+7]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+6] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+6], lfr_bp2[i*NB_BYTES_BP2+6]);
             #endif
                     // S55
                     significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent);  // 0.5 <= significand < 1
                                                                                                                // S55 = significand * 2^exponent
             #ifdef DEBUG_TCH
                     printf("S55         : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);
                     printf("significand : %16.8e\n",significand);
                     printf("exponent    : %d\n"    ,exponent);
             #endif
                     if (exponent < expmin) { // value should be >= 0.5 * 2^expmin
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     if (exponent > expmax) { // value should be <  0.5 * 2^(expmax+1)
                       exponent = expmax;
                       significand = 1.0;     // max value that can be recorded
                     }
                     if (significand == 0) {  // in that case exponent == 0 too
                       exponent = expmin;
                       significand = 0.5;     // min value that can be recorded
                     }
                     autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding
                                                                              // where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)
                     exp = (uint16_t) (exponent-expmin);      // Shift and cast into a 16-bit unsigned int where just
                                                              // the first nbitexp bits are used (0, ..., 2^nbitexp-1)
                     tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the
                                                              // left place of the significand bits (nbitsig),
                                                              // making the 16-bit word to be recorded
                     pt_uint8 = (uint8_t*) &tmp_uint16;       // Affect an uint8_t pointer with the adress of tmp_uint16
             #ifdef MSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record LSB of tmp_uint16
                     //printf("MSB:\n");
             #endif
             #ifdef LSB_FIRST_TCH
                     lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record MSB of tmp_uint16
                     lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record LSB of tmp_uint16
                     //printf("LSB:\n");
             #endif
             #ifdef DEBUG_TCH
                     printf("autocor for S55 significand : %u\n",autocor);
                     printf("exp for S55 exponent : %u\n",exp);
                     printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);
                     printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+9] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+9], lfr_bp2[i*NB_BYTES_BP2+9]);
                     printf("lfr_bp2[i*NB_BYTES_BP2+8] : %3u or %2x\n",lfr_bp2[i*NB_BYTES_BP2+8], lfr_bp2[i*NB_BYTES_BP2+8]);
             #endif
                 }
             }

             // In the frame of RPW LFR Sofware ICD Issue1 Rev8 (05/07/2013)
             // version 1.0: 31/07/2013
             // version 1.1: 02/04/2014
             // version 1.2: 30/04/2014
             // version 1.3: 02/05/2014
             // version 1.4: 16/05/2014
             // version 1.5: 20/05/2014
             // version 1.6: 19/12/2014
             // version 1.7: 15/01/2015 (modifs de Paul + correction erreurs qui se compensaient (LSB <=> MSB + indices [0,2] <=> [1,3])
             #ifndef BASIC_PARAMETERS_H_INCLUDED
             #define BASIC_PARAMETERS_H_INCLUDED
             #include <math.h>
             #include <stdio.h>
             #include <stdint.h>
             #include "basic_parameters_params.h"
             void BP1_set(float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1);
+            float sqrt_ple( float x );
             void BP2_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2);
+            void BP3_set(float * compressed_spec_mat, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp2);
             void init_k_coefficients( float *k_coefficients_f0, unsigned char nb_binscompressed_matrix );
             //***********************************
             // STATIC INLINE FUNCTION DEFINITIONS
             #endif // BASIC_PARAMETERS_H_INCLUDED