HG_REPOSITORIES/LPP/INSTRUMENTATION/USERS/CHUST/LFR_basic-parameters Commit - r16:a586fe639ac1

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

// version 1.8: 02/02/2015 (gestion des divisions par zéro)

10

11

#ifndef BASIC_PARAMETERS_H_INCLUDED

12

#ifndef BASIC_PARAMETERS_H_INCLUDED

12

#define BASIC_PARAMETERS_H_INCLUDED

13

#define BASIC_PARAMETERS_H_INCLUDED

13

14

#include <math.h>

15

#include <math.h>

15

#include <stdio.h>

16

#include <stdio.h>

16

#include <stdint.h>

17

#include <stdint.h>

17

18

#include "basic_parameters_params.h"

19

#include "basic_parameters_params.h"

19

20

static inline void BP1_set(float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1);

21

static inline void BP1_set(float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1);

21

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

22

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

22

23

void init_k_coefficients( float *k_coefficients_f0, unsigned char nb_binscompressed_matrix );

24

void init_k_coefficients( float *k_coefficients_f0, unsigned char nb_binscompressed_matrix );

24

25

//***********************************

26

//***********************************

26

// STATIC INLINE FUNCTION DEFINITIONS

27

// STATIC INLINE FUNCTION DEFINITIONS

27

28

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

29

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

29

float PSDB; // 32-bit floating point

30

float PSDB; // 32-bit floating point

30

float PSDE;

31

float PSDE;

31

float tmp;

32

float tmp;

32

float NVEC_V0;

33

float NVEC_V0;

33

float NVEC_V1;

34

float NVEC_V1;

34

float NVEC_V2;

35

float NVEC_V2;

35

float aux;

36

float aux;

36

float tr_SB_SB;

37

float tr_SB_SB;

37

float e_cross_b_re;

38

float e_cross_b_re;

38

float e_cross_b_im;

39

float e_cross_b_im;

39

float n_cross_e_scal_b_re;

40

float n_cross_e_scal_b_re;

40

float n_cross_e_scal_b_im;

41

float n_cross_e_scal_b_im;

41

float ny;

42

float ny;

42

float nz;

43

float nz;

43

float bx_bx_star;

44

float bx_bx_star;

44

float vphi;

45

float vphi;

45

float significand;

46

float significand;

46

int exponent; // 32-bit signed integer

47

int exponent; // 32-bit signed integer

47

float alpha_M;

48

float alpha_M;

48

49

uint8_t nbitexp; // 8-bit unsigned integer

50

uint8_t nbitexp; // 8-bit unsigned integer

50

uint8_t nbitsig;

51

uint8_t nbitsig;

51

uint8_t tmp_uint8;

52

uint8_t tmp_uint8;

52

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

53

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

53

int8_t expmin; // 8-bit signed integer

54

int8_t expmin; // 8-bit signed integer

54

int8_t expmax;

55

int8_t expmax;

55

uint16_t rangesig; // 16-bit unsigned integer

56

uint16_t rangesig; // 16-bit unsigned integer

56

uint16_t psd;

57

uint16_t psd;

57

uint16_t exp;

58

uint16_t exp;

58

uint16_t tmp_uint16;

59

uint16_t tmp_uint16;

59

uint16_t i;

60

uint16_t i;

60

61

alpha_M = 45 * (3.1415927/180);

62

alpha_M = 45 * (3.1415927/180);

62

63

#ifdef DEBUG_TCH

64

#ifdef DEBUG_TCH

64

printf("BP1 : \n");

65

printf("BP1 : \n");

65

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

66

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

66

#endif

67

#endif

67

68

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

69

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

69

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

70

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

70

expmax = 30; // maximum value of the exponent

71

expmax = 30; // maximum value of the exponent

71

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

72

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

72

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

73

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

73

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

74

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

74

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

75

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

75

76

#ifdef DEBUG_TCH

77

#ifdef DEBUG_TCH

77

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

78

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

78

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

79

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

79

#endif

80

#endif

80

81

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

82

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

82

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

83

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

83

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

84

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

84

PSDB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] // S11

85

PSDB = compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX] // S11

85

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

86

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

86

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

87

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

87

88

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

89

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

89

// PSDB/3 = significand * 2^exponent

90

// PSDB/3 = significand * 2^exponent

90

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

91

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

91

92

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

93

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

93

exponent = expmin;

94

exponent = expmin;

94

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

95

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

95

}

96

}

96

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

97

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

97

exponent = expmax;

98

exponent = expmax;

98

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

99

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

99

}

100

}

100

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

101

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

101

exponent = expmin;

102

exponent = expmin;

102

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

103

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

103

}

104

}

104

105

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

106

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

106

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

107

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

107

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

108

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

108

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

109

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

109

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

110

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

110

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

111

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

111

// making the 16-bit word to be recorded

112

// making the 16-bit word to be recorded

112

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

113

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

113

#ifdef MSB_FIRST_TCH

114

#ifdef MSB_FIRST_TCH

114

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

115

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

115

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

116

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

116

#endif

117

#endif

117

#ifdef LSB_FIRST_TCH

118

#ifdef LSB_FIRST_TCH

118

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

119

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

119

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

120

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

120

#endif

121

#endif

121

#ifdef DEBUG_TCH

122

#ifdef DEBUG_TCH

122

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

123

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

123

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

124

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

124

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

125

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

125

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

126

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

126

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

127

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

127

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

128

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

128

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

129

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

129

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

130

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

130

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]);

131

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]);

131

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]);

132

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]);

132

#endif

133

#endif

133

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

134

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

134

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

135

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

135

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

136

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

136

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

137

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

137

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

138

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

138

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

139

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

139

140

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

141

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

141

// PSDE/2 = significand * 2^exponent

142

// PSDE/2 = significand * 2^exponent

142

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

143

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

143

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

144

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

144

145

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

146

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

146

exponent = expmin;

147

exponent = expmin;

147

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

148

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

148

}

149

}

149

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

150

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

150

exponent = expmax;

151

exponent = expmax;

151

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

152

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

152

}

153

}

153

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

154

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

154

exponent = expmin;

155

exponent = expmin;

155

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

156

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

156

}

157

}

157

158

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

159

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

159

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

160

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

160

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

161

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

161

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

162

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

162

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

163

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

163

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

164

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

164

// making the 16-bit word to be recorded

165

// making the 16-bit word to be recorded

165

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

166

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

166

#ifdef MSB_FIRST_TCH

167

#ifdef MSB_FIRST_TCH

167

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

168

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

168

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

169

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

169

#endif

170

#endif

170

#ifdef LSB_FIRST_TCH

171

#ifdef LSB_FIRST_TCH

171

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

172

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

172

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

173

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

173

#endif

174

#endif

174

#ifdef DEBUG_TCH

175

#ifdef DEBUG_TCH

175

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

176

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

176

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

177

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

177

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

178

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

178

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

179

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

179

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

180

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

180

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

181

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

181

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

182

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

182

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

183

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

183

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]);

184

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]);

184

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]);

185

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]);

185

#endif

186

#endif

186

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

187

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

187

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

188

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

188

// == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits

189

// == PA_LFR_SC_BP1_NVEC_V1_F0 == 8 bits

189

// == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit

190

// == PA_LFR_SC_BP1_NVEC_V2_F0 == 1 sign bit

190

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

191

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

191

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

192

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

192

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

193

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

193

);

194

);

195

if (tmp != 0.) { // no division by 0.

194

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

196

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

195

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

197

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

196

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

198

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

197

199

}

200

else

201

{

202

NVEC_V0 = 0.;

203

NVEC_V1 = 0.;

204

NVEC_V2 = 0.;

205

}

198

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

206

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

199

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

207

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

200

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

208

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

201

#ifdef LSB_FIRST_TCH

209

#ifdef LSB_FIRST_TCH

202

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)

210

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)

203

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

211

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

204

#endif

212

#endif

205

#ifdef MSB_FIRST_TCH

213

#ifdef MSB_FIRST_TCH

206

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)

214

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)

207

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

215

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

208

#endif

216

#endif

209

#ifdef DEBUG_TCH

217

#ifdef DEBUG_TCH

210

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

218

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

211

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

219

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

212

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

220

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

213

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

221

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

214

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

222

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

215

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

223

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

216

#endif

224

#endif

217

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

225

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

218

// BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits

226

// BP1 ellipticity == PA_LFR_SC_BP1_ELLIP_F0 == 4 bits

227

if (PSDB != 0.) { // no division by 0.

219

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

228

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

220

229

}

230

else

231

{

232

aux = 0.;

233

}

221

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

234

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

222

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

235

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

223

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

236

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

224

// of the sign bit of NVEC_V2 (recorded

237

// of the sign bit of NVEC_V2 (recorded

225

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

238

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

226

#ifdef DEBUG_TCH

239

#ifdef DEBUG_TCH

227

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

240

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

228

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

241

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

229

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

242

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

230

#endif

243

#endif

231

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

244

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

232

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

245

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

233

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

246

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

234

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

247

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

235

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

248

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

236

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

249

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

237

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

250

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

238

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

251

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

239

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

252

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

240

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

253

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

241

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

254

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

242

aux = PSDB*PSDB;

255

aux = PSDB*PSDB;

256

if (aux != 0.) { // no division by 0.

243

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

257

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

244

258

}

259

else

260

{

261

tmp = 0.;

262

}

245

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

263

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

246

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

264

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

247

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

265

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

248

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

266

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

249

#ifdef DEBUG_TCH

267

#ifdef DEBUG_TCH

250

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

268

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

251

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

269

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

252

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

270

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

253

#endif

271

#endif

254

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

272

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

255

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

273

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

256

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

274

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

257

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

275

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

258

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

276

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

259

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

277

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

260

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

278

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

261

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

279

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

262

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

280

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

263

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

281

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

264

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

282

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

265

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

283

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

266

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

284

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

267

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

285

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

268

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

286

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

269

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

287

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

270

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

288

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

271

// k_ji = k_ij

289

// k_ji = k_ij

272

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

290

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

273

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

291

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

274

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

292

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

275

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

293

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

276

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

294

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

277

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

295

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

278

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

296

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

279

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

297

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

280

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

298

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

281

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

299

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

282

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

300

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

283

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

301

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

284

#ifdef DEBUG_TCH

302

#ifdef DEBUG_TCH

285

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

303

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

286

#endif

304

#endif

287

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

305

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

288

#ifdef LSB_FIRST_TCH

306

#ifdef LSB_FIRST_TCH

289

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)

307

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)

290

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

308

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

291

// of lfr_bp1[i*NB_BYTES_BP1+0]

309

// of lfr_bp1[i*NB_BYTES_BP1+0]

292

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

310

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

293

#endif

311

#endif

294

#ifdef MSB_FIRST_TCH

312

#ifdef MSB_FIRST_TCH

295

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)

313

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)

296

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

314

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

297

// of lfr_bp1[i*NB_BYTES_BP1+0]

315

// of lfr_bp1[i*NB_BYTES_BP1+0]

298

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

316

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

299

#endif

317

#endif

300

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

318

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

301

// ReaSX/2 = significand * 2^exponent

319

// ReaSX/2 = significand * 2^exponent

302

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

320

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

303

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

321

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

304

322

305

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

323

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

306

exponent = expmin;

324

exponent = expmin;

307

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

325

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

308

}

326

}

309

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

327

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

310

exponent = expmax;

328

exponent = expmax;

311

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

329

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

312

}

330

}

313

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

331

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

314

exponent = expmin;

332

exponent = expmin;

315

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

333

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

316

}

334

}

317

335

318

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

336

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

319

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

337

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

320

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

338

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

321

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

339

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

322

#ifdef DEBUG_TCH

340

#ifdef DEBUG_TCH

323

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

341

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

324

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

342

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

325

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

343

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

326

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

344

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

327

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

345

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

328

#endif

346

#endif

329

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

347

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

330

// with lfr_bp1[i*NB_BYTES_BP1+7]

348

// with lfr_bp1[i*NB_BYTES_BP1+7]

331

#ifdef DEBUG_TCH

349

#ifdef DEBUG_TCH

332

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

350

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

333

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

351

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

334

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

352

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

335

#endif

353

#endif

336

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

354

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

337

#ifdef LSB_FIRST_TCH

355

#ifdef LSB_FIRST_TCH

338

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)

356

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)

339

#endif

357

#endif

340

#ifdef MSB_FIRST_TCH

358

#ifdef MSB_FIRST_TCH

341

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)

359

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)

342

#endif

360

#endif

343

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

361

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

344

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

362

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

345

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

363

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

346

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

364

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

347

#ifdef DEBUG_TCH

365

#ifdef DEBUG_TCH

348

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

366

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

349

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

367

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

350

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

368

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

351

#endif

369

#endif

352

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

370

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

353

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

371

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

354

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

372

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

355

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

373

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

356

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

374

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

357

nz = NVEC_V0;

375

nz = NVEC_V0;

358

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

376

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

359

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

377

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

360

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

378

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

361

379

362

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

380

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

363

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

381

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

364

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

382

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

365

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

383

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

366

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

384

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

367

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

385

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

368

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

386

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

369

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

387

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

370

+ 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

388

+ 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

371

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

389

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

372

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

390

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

373

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

391

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

374

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

392

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

375

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

393

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

376

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

394

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

377

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

395

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

378

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

396

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

379

// k_ji = k_ij

397

// k_ji = k_ij

380

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

398

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

381

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

399

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

382

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

400

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

383

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

401

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

384

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

402

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

385

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

403

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

386

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

404

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

387

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

405

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

388

+ 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

406

+ 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

389

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

407

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

390

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

408

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

391

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

409

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

392

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

410

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

393

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

411

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

394

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

412

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

395

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

413

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

396

#ifdef DEBUG_TCH

414

#ifdef DEBUG_TCH

397

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

415

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

398

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

416

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

399

#endif

417

#endif

400

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

418

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

401

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

419

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

402

#ifdef LSB_FIRST_TCH

420

#ifdef LSB_FIRST_TCH

403

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)

421

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)

404

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

422

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

405

// of lfr_bp1[i*NB_BYTES_BP1+2]

423

// of lfr_bp1[i*NB_BYTES_BP1+2]

406

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|

424

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|

407

#endif

425

#endif

408

#ifdef MSB_FIRST_TCH

426

#ifdef MSB_FIRST_TCH

409

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)

427

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)

410

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

428

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

411

// of lfr_bp1[i*NB_BYTES_BP1+2]

429

// of lfr_bp1[i*NB_BYTES_BP1+2]

412

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|

430

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|

413

#endif

431

#endif

432

if (bx_bx_star != 0.) { // no division by 0.

414

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

433

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

415

434

}

435

else

436

{

437

vphi = 1.e+20; // Put a huge value

438

}

416

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

439

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

417

// vphi/2 = significand * 2^exponent

440

// vphi/2 = significand * 2^exponent

418

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

441

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

419

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

442

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

420

443

421

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

444

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

422

exponent = expmin;

445

exponent = expmin;

423

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

446

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

424

}

447

}

425

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

448

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

426

exponent = expmax;

449

exponent = expmax;

427

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

450

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

428

}

451

}

429

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

452

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

430

exponent = expmin;

453

exponent = expmin;

431

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

454

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

432

}

455

}

433

#ifdef DEBUG_TCH

456

#ifdef DEBUG_TCH

434

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

457

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

435

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

458

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

436

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

459

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

437

#endif

460

#endif

438

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

461

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

439

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

462

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

440

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

463

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

441

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

464

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

442

#ifdef DEBUG_TCH

465

#ifdef DEBUG_TCH

443

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

466

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

444

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

467

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

445

#endif

468

#endif

446

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

469

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

447

// with lfr_bp1[i*NB_BYTES_BP1+8]

470

// with lfr_bp1[i*NB_BYTES_BP1+8]

448

#ifdef DEBUG_TCH

471

#ifdef DEBUG_TCH

449

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

472

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

450

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

473

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

451

#endif

474

#endif

452

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

475

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

453

#ifdef LSB_FIRST_TCH

476

#ifdef LSB_FIRST_TCH

454

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)

477

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)

455

#endif

478

#endif

456

#ifdef MSB_FIRST_TCH

479

#ifdef MSB_FIRST_TCH

457

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)

480

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)

458

#endif

481

#endif

459

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

482

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

460

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

483

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

461

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

484

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

462

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

485

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

463

#ifdef DEBUG_TCH

486

#ifdef DEBUG_TCH

464

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

487

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

465

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

488

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

466

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

489

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

467

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

490

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

468

#endif

491

#endif

469

}

492

}

470

}

493

}

471

494

472

void BP2_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * lfr_bp2 )

495

void BP2_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * lfr_bp2 )

473

{

496

{

474

float cross_re; // 32-bit floating point

497

float cross_re; // 32-bit floating point

475

float cross_im;

498

float cross_im;

476

float aux;

499

float aux;

477

float significand;

500

float significand;

478

int exponent; // 32-bit signed integer

501

int exponent; // 32-bit signed integer

479

uint8_t nbitexp; // 8-bit unsigned integer

502

uint8_t nbitexp; // 8-bit unsigned integer

480

uint8_t nbitsig;

503

uint8_t nbitsig;

481

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

504

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

482

int8_t expmin; // 8-bit signed integer

505

int8_t expmin; // 8-bit signed integer

483

int8_t expmax;

506

int8_t expmax;

484

uint16_t rangesig; // 16-bit unsigned integer

507

uint16_t rangesig; // 16-bit unsigned integer

485

uint16_t autocor;

508

uint16_t autocor;

486

uint16_t exp;

509

uint16_t exp;

487

uint16_t tmp_uint16;

510

uint16_t tmp_uint16;

488

uint16_t i;

511

uint16_t i;

489

512

490

#ifdef DEBUG_TCH

513

#ifdef DEBUG_TCH

491

printf("BP2 : \n");

514

printf("BP2 : \n");

492

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

515

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

493

#endif

516

#endif

494

517

495

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

518

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

496

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

519

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

497

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

520

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

498

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

521

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

499

expmax = 32;

522

expmax = 32;

500

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

523

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

501

524

502

#ifdef DEBUG_TCH

525

#ifdef DEBUG_TCH

503

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

526

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

504

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

527

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

505

#endif

528

#endif

506

529

507

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

530

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

508

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

531

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

509

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

532

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

510

// == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits

533

// == PA_LFR_SC_BP2_CROSS_RE_0_F0 == 8 bits

511

// == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits

534

// == PA_LFR_SC_BP2_CROSS_IM_0_F0 == 8 bits

512

// == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits

535

// == PA_LFR_SC_BP2_CROSS_RE_1_F0 == 8 bits

513

// == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits

536

// == PA_LFR_SC_BP2_CROSS_IM_1_F0 == 8 bits

514

// == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits

537

// == PA_LFR_SC_BP2_CROSS_RE_2_F0 == 8 bits

515

// == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits

538

// == PA_LFR_SC_BP2_CROSS_IM_2_F0 == 8 bits

516

// == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits

539

// == PA_LFR_SC_BP2_CROSS_RE_3_F0 == 8 bits

517

// == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits

540

// == PA_LFR_SC_BP2_CROSS_IM_3_F0 == 8 bits

518

// == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits

541

// == PA_LFR_SC_BP2_CROSS_RE_4_F0 == 8 bits

519

// == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits

542

// == PA_LFR_SC_BP2_CROSS_IM_4_F0 == 8 bits

520

// == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits

543

// == PA_LFR_SC_BP2_CROSS_RE_5_F0 == 8 bits

521

// == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits

544

// == PA_LFR_SC_BP2_CROSS_IM_5_F0 == 8 bits

522

// == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits

545

// == PA_LFR_SC_BP2_CROSS_RE_6_F0 == 8 bits

523

// == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits

546

// == PA_LFR_SC_BP2_CROSS_IM_6_F0 == 8 bits

524

// == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits

547

// == PA_LFR_SC_BP2_CROSS_RE_7_F0 == 8 bits

525

// == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits

548

// == PA_LFR_SC_BP2_CROSS_IM_7_F0 == 8 bits

526

// == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits

549

// == PA_LFR_SC_BP2_CROSS_RE_8_F0 == 8 bits

527

// == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits

550

// == PA_LFR_SC_BP2_CROSS_IM_8_F0 == 8 bits

528

// == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits

551

// == PA_LFR_SC_BP2_CROSS_RE_9_F0 == 8 bits

529

// == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits

552

// == PA_LFR_SC_BP2_CROSS_IM_9_F0 == 8 bits

530

// S12

553

// S12

531

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

554

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

555

if (aux != 0.) { // no division by 0.

532

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

556

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

533

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

557

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

558

}

559

else

560

{

561

cross_re = 0.;

562

cross_im = 0.;

563

}

534

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

564

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

535

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

565

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

536

#ifdef DEBUG_TCH

566

#ifdef DEBUG_TCH

537

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

567

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

538

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

568

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

539

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

569

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

540

#endif

570

#endif

541

// S13

571

// S13

542

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

572

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

573

if (aux != 0.) { // no division by 0.

543

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

574

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

544

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

575

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

576

}

577

else

578

{

579

cross_re = 0.;

580

cross_im = 0.;

581

}

545

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

582

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

546

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

583

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

547

#ifdef DEBUG_TCH

584

#ifdef DEBUG_TCH

548

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

585

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

549

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

586

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

550

#endif

587

#endif

551

// S14

588

// S14

552

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

589

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

590

if (aux != 0.) { // no division by 0.

553

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

591

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

554

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

592

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

593

}

594

else

595

{

596

cross_re = 0.;

597

cross_im = 0.;

598

}

555

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

599

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

556

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

600

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

557

#ifdef DEBUG_TCH

601

#ifdef DEBUG_TCH

558

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

602

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

559

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

603

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

560

#endif

604

#endif

561

// S15

605

// S15

562

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

606

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

607

if (aux != 0.) { // no division by 0.

563

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

608

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

564

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

609

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

610

}

611

else

612

{

613

cross_re = 0.;

614

cross_im = 0.;

615

}

565

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

566

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

567

#ifdef DEBUG_TCH

618

#ifdef DEBUG_TCH

568

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]);

569

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]);

570

#endif

621

#endif

571

// S23

622

// S23

572

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(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+9]*compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+16]);

624

if (aux != 0.) { // no division by 0.

573

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

625

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

574

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

626

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

627

}

628

else

629

{

630

cross_re = 0.;

631

cross_im = 0.;

632

}

575

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

633

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

576

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

634

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

577

#ifdef DEBUG_TCH

635

#ifdef DEBUG_TCH

578

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

636

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

579

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

637

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

580

#endif

638

#endif

581

// S24

639

// S24

582

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

640

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

641

if (aux != 0.) { // no division by 0.

583

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

642

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

584

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

643

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

644

}

645

else

646

{

647

cross_re = 0.;

648

cross_im = 0.;

649

}

585

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

650

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

586

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

651

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

587

#ifdef DEBUG_TCH

652

#ifdef DEBUG_TCH

588

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

653

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

589

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

654

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

590

#endif

655

#endif

591

// S25

656

// S25

592

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

657

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

658

if (aux != 0.) { // no division by 0.

593

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

659

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

594

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

660

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

661

}

662

else

663

{

664

cross_re = 0.;

665

cross_im = 0.;

666

}

595

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

667

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

596

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

668

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

597

#ifdef DEBUG_TCH

669

#ifdef DEBUG_TCH

598

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

670

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

599

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

671

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

600

#endif

672

#endif

601

// S34

673

// S34

602

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

674

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

675

if (aux != 0.) { // no division by 0.

603

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

676

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

604

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

677

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

678

}

679

else

680

{

681

cross_re = 0.;

682

cross_im = 0.;

683

}

605

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

684

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

606

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

685

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

607

#ifdef DEBUG_TCH

686

#ifdef DEBUG_TCH

608

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

687

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

609

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

688

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

610

#endif

689

#endif

611

// S35

690

// S35

612

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

691

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

692

if (aux != 0.) { // no division by 0.

613

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

693

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

614

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

694

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

695

}

696

else

697

{

698

cross_re = 0.;

699

cross_im = 0.;

700

}

615

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

701

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

616

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

702

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

617

#ifdef DEBUG_TCH

703

#ifdef DEBUG_TCH

618

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

704

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

619

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

705

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

620

#endif

706

#endif

621

// S45

707

// S45

622

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

708

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

709

if (aux != 0.) { // no division by 0.

623

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

710

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

624

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

711

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

712

}

713

else

714

{

715

cross_re = 0.;

716

cross_im = 0.;

717

}

625

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

718

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

626

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

719

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

627

#ifdef DEBUG_TCH

720

#ifdef DEBUG_TCH

628

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

721

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

629

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

722

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

630

#endif

723

#endif

631

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

724

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

632

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

725

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

633

// == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits

726

// == PA_LFR_SC_BP2_AUTO_A0_F0 == 16 bits

634

// == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits

727

// == PA_LFR_SC_BP2_AUTO_A1_F0 == 16 bits

635

// == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits

728

// == PA_LFR_SC_BP2_AUTO_A2_F0 == 16 bits

636

// == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits

729

// == PA_LFR_SC_BP2_AUTO_A3_F0 == 16 bits

637

// == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits

730

// == PA_LFR_SC_BP2_AUTO_A4_F0 == 16 bits

638

// S11

731

// S11

639

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

732

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

640

// S11 = significand * 2^exponent

733

// S11 = significand * 2^exponent

641

#ifdef DEBUG_TCH

734

#ifdef DEBUG_TCH

642

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

735

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

643

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

736

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

644

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

737

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

645

#endif

738

#endif

646

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

739

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

647

exponent = expmin;

740

exponent = expmin;

648

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

741

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

649

}

742

}

650

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

743

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

651

exponent = expmax;

744

exponent = expmax;

652

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

745

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

653

}

746

}

654

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

747

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

655

exponent = expmin;

748

exponent = expmin;

656

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

749

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

657

}

750

}

658

751

659

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

752

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

660

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

753

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

661

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

754

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

662

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

755

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

663

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

756

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

664

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

757

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

665

// making the 16-bit word to be recorded

758

// making the 16-bit word to be recorded

666

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

759

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

667

#ifdef MSB_FIRST_TCH

760

#ifdef MSB_FIRST_TCH

668

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

761

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

669

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

762

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

670

#endif

763

#endif

671

#ifdef LSB_FIRST_TCH

764

#ifdef LSB_FIRST_TCH

672

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

765

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

673

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

766

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

674

#endif

767

#endif

675

#ifdef DEBUG_TCH

768

#ifdef DEBUG_TCH

676

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

769

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

677

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

770

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

678

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

771

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

679

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

772

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

680

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]);

773

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]);

681

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]);

774

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]);

682

#endif

775

#endif

683

// S22

776

// S22

684

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

777

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

685

// S22 = significand * 2^exponent

778

// S22 = significand * 2^exponent

686

#ifdef DEBUG_TCH

779

#ifdef DEBUG_TCH

687

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

780

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

688

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

781

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

689

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

782

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

690

#endif

783

#endif

691

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

784

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

692

exponent = expmin;

785

exponent = expmin;

693

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

786

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

694

}

787

}

695

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

788

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

696

exponent = expmax;

789

exponent = expmax;

697

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

790

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

698

}

791

}

699

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

792

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

700

exponent = expmin;

793

exponent = expmin;

701

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

794

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

702

}

795

}

703

796

704

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

797

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

705

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

798

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

706

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

799

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

707

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

800

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

708

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

801

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

709

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

802

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

710

// making the 16-bit word to be recorded

803

// making the 16-bit word to be recorded

711

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

804

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

712

#ifdef MSB_FIRST_TCH

805

#ifdef MSB_FIRST_TCH

713

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

806

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

714

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

807

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

715

#endif

808

#endif

716

#ifdef LSB_FIRST_TCH

809

#ifdef LSB_FIRST_TCH

717

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

810

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

718

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

811

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

719

#endif

812

#endif

720

#ifdef DEBUG_TCH

813

#ifdef DEBUG_TCH

721

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

814

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

722

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

815

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

723

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

816

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

724

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

817

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

725

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]);

818

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]);

726

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]);

819

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]);

727

#endif

820

#endif

728

// S33

821

// S33

729

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

822

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

730

// S33 = significand * 2^exponent

823

// S33 = significand * 2^exponent

731

#ifdef DEBUG_TCH

824

#ifdef DEBUG_TCH

732

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

825

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

733

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

826

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

734

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

827

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

735

#endif

828

#endif

736

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

829

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

737

exponent = expmin;

830

exponent = expmin;

738

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

831

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

739

}

832

}

740

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

833

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

741

exponent = expmax;

834

exponent = expmax;

742

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

835

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

743

}

836

}

744

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

837

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

745

exponent = expmin;

838

exponent = expmin;

746

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

839

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

747

}

840

}

748

841

749

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

842

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

750

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

843

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

751

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

844

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

752

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

845

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

753

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

846

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

754

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

847

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

755

// making the 16-bit word to be recorded

848

// making the 16-bit word to be recorded

756

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

849

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

757

#ifdef MSB_FIRST_TCH

850

#ifdef MSB_FIRST_TCH

758

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

851

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

759

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

852

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

760

#endif

853

#endif

761

#ifdef LSB_FIRST_TCH

854

#ifdef LSB_FIRST_TCH

762

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

855

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

763

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

856

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

764

#endif

857

#endif

765

#ifdef DEBUG_TCH

858

#ifdef DEBUG_TCH

766

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

859

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

767

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

860

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

768

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

861

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

769

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

862

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

770

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]);

863

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]);

771

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]);

864

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]);

772

#endif

865

#endif

773

// S44

866

// S44

774

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

867

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

775

// S44 = significand * 2^exponent

868

// S44 = significand * 2^exponent

776

#ifdef DEBUG_TCH

869

#ifdef DEBUG_TCH

777

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

870

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

778

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

871

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

779

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

872

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

780

#endif

873

#endif

781

874

782

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

875

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

783

exponent = expmin;

876

exponent = expmin;

784

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

877

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

785

}

878

}

786

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

879

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

787

exponent = expmax;

880

exponent = expmax;

788

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

881

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

789

}

882

}

790

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

883

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

791

exponent = expmin;

884

exponent = expmin;

792

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

885

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

793

}

886

}

794

887

795

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

888

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

796

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

889

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

797

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

890

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

798

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

891

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

799

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

892

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

800

// left place of the significand bits (nbitsig),

893

// left place of the significand bits (nbitsig),

801

// making the 16-bit word to be recorded

894

// making the 16-bit word to be recorded

802

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

895

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

803

#ifdef MSB_FIRST_TCH

896

#ifdef MSB_FIRST_TCH

804

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record MSB of tmp_uint16

897

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[0]; // Record MSB of tmp_uint16

805

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record LSB of tmp_uint16

898

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[1]; // Record LSB of tmp_uint16

806

#endif

899

#endif

807

#ifdef LSB_FIRST_TCH

900

#ifdef LSB_FIRST_TCH

808

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record MSB of tmp_uint16

901

lfr_bp2[i*NB_BYTES_BP2+6] = pt_uint8[1]; // Record MSB of tmp_uint16

809

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record LSB of tmp_uint16

902

lfr_bp2[i*NB_BYTES_BP2+7] = pt_uint8[0]; // Record LSB of tmp_uint16

810

#endif

903

#endif

811

#ifdef DEBUG_TCH

904

#ifdef DEBUG_TCH

812

printf("autocor for S44 significand : %u\n",autocor);

905

printf("autocor for S44 significand : %u\n",autocor);

813

printf("exp for S44 exponent : %u\n",exp);

906

printf("exp for S44 exponent : %u\n",exp);

814

printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

907

printf("pt_uint8[1] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

815

printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

908

printf("pt_uint8[0] for S44 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

816

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]);

909

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]);

817

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]);

910

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]);

818

#endif

911

#endif

819

// S55

912

// S55

820

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent); // 0.5 <= significand < 1

913

significand = frexpf(compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24], &exponent); // 0.5 <= significand < 1

821

// S55 = significand * 2^exponent

914

// S55 = significand * 2^exponent

822

#ifdef DEBUG_TCH

915

#ifdef DEBUG_TCH

823

printf("S55 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

916

printf("S55 : %16.8e\n",compressed_spec_mat[i*NB_VALUES_PER_SPECTRAL_MATRIX+24]);

824

printf("significand : %16.8e\n",significand);

917

printf("significand : %16.8e\n",significand);

825

printf("exponent : %d\n" ,exponent);

918

printf("exponent : %d\n" ,exponent);

826

#endif

919

#endif

827

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

920

if (exponent < expmin) { // value should be >= 0.5 * 2^expmin

828

exponent = expmin;

921

exponent = expmin;

829

significand = 0.5; // min value that can be recorded

922

significand = 0.5; // min value that can be recorded

830

}

923

}

831

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

924

if (exponent > expmax) { // value should be < 0.5 * 2^(expmax+1)

832

exponent = expmax;

925

exponent = expmax;

833

significand = 1.0; // max value that can be recorded

926

significand = 1.0; // max value that can be recorded

834

}

927

}

835

if (significand == 0) { // in that case exponent == 0 too

928

if (significand == 0) { // in that case exponent == 0 too

836

exponent = expmin;

929

exponent = expmin;

837

significand = 0.5; // min value that can be recorded

930

significand = 0.5; // min value that can be recorded

838

}

931

}

839

932

840

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

933

autocor = (uint16_t) ((significand*2-1)*rangesig + 0.5); // Shift and cast into a 16-bit unsigned int with rounding

841

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

934

// where just the first nbitsig bits are used (0, ..., 2^nbitsig-1)

842

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

935

exp = (uint16_t) (exponent-expmin); // Shift and cast into a 16-bit unsigned int where just

843

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

936

// the first nbitexp bits are used (0, ..., 2^nbitexp-1)

844

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

937

tmp_uint16 = autocor | (exp << nbitsig); // Put the exponent bits (nbitexp) next to the

845

// left place of the significand bits (nbitsig),

938

// left place of the significand bits (nbitsig),

846

// making the 16-bit word to be recorded

939

// making the 16-bit word to be recorded

847

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

940

pt_uint8 = (uint8_t*) &tmp_uint16; // Affect an uint8_t pointer with the adress of tmp_uint16

848

#ifdef MSB_FIRST_TCH

941

#ifdef MSB_FIRST_TCH

849

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record MSB of tmp_uint16

942

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[0]; // Record MSB of tmp_uint16

850

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record LSB of tmp_uint16

943

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[1]; // Record LSB of tmp_uint16

851

//printf("MSB:\n");

944

//printf("MSB:\n");

852

#endif

945

#endif

853

#ifdef LSB_FIRST_TCH

946

#ifdef LSB_FIRST_TCH

854

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record MSB of tmp_uint16

947

lfr_bp2[i*NB_BYTES_BP2+8] = pt_uint8[1]; // Record MSB of tmp_uint16

855

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record LSB of tmp_uint16

948

lfr_bp2[i*NB_BYTES_BP2+9] = pt_uint8[0]; // Record LSB of tmp_uint16

856

//printf("LSB:\n");

949

//printf("LSB:\n");

857

#endif

950

#endif

858

#ifdef DEBUG_TCH

951

#ifdef DEBUG_TCH

859

printf("autocor for S55 significand : %u\n",autocor);

952

printf("autocor for S55 significand : %u\n",autocor);

860

printf("exp for S55 exponent : %u\n",exp);

953

printf("exp for S55 exponent : %u\n",exp);

861

printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

954

printf("pt_uint8[1] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[1], pt_uint8[1]);

862

printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

955

printf("pt_uint8[0] for S55 exponent + significand : %.3d or %2x\n",pt_uint8[0], pt_uint8[0]);

863

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]);

956

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]);

864

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]);

957

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]);

865

#endif

958

#endif

866

}

959

}

867

}

960

}

868

961

869

962

870

#endif // BASIC_PARAMETERS_H_INCLUDED

963

#endif // BASIC_PARAMETERS_H_INCLUDED

	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])
+            // version 1.8: 02/02/2015 (gestion des divisions par zéro)
             #ifndef BASIC_PARAMETERS_H_INCLUDED
             #define BASIC_PARAMETERS_H_INCLUDED
             #include <math.h>
             #include <stdio.h>
             #include <stdint.h>
             #include "basic_parameters_params.h"
             static inline void BP1_set(float * compressed_spec_mat, float * k_coeff_intercalib, unsigned char nb_bins_compressed_spec_mat, unsigned char * lfr_bp1);
             static inline void BP2_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
             void BP1_set( float * compressed_spec_mat, float * k_coeff_intercalib, uint8_t nb_bins_compressed_spec_mat, uint8_t * 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
                                +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
                                );
+                    if (tmp != 0.) { // no division by 0.
                         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
+                    }
+                    else
+                    {
+                        NVEC_V0 = 0.;
+                        NVEC_V1 = 0.;
+                        NVEC_V2 = 0.;
+                    }
                     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
+                    if (PSDB != 0.) { // no division by 0.
                         aux = 2*tmp / PSDB;                   // Compute the ellipticity
+                    }
+                    else
+                    {
+                        aux = 0.;
+                    }
                     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;
+                    if (aux != 0.) { // no division by 0.
                         tmp = ( 3*tr_SB_SB - aux ) / ( 2 * aux );  // Compute the degree of polarisation
+                    }
+                    else
+                    {
+                        tmp = 0.;
+                    }
                     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
+                    if (bx_bx_star != 0.) { // no division by 0.
                         vphi = n_cross_e_scal_b_re / bx_bx_star;  // Compute |VPHI|
+                    }
+                    else
+                    {
+                        vphi = 1.e+20;                         // Put a huge value
+                    }
                     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
                 }
             }
             void BP2_set( float * compressed_spec_mat, uint8_t nb_bins_compressed_spec_mat, uint8_t * 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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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]);
+                    if (aux != 0.) { // no division by 0.
                     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;
+                    }
+                    else
+                    {
+                    cross_re = 0.;
+                    cross_im = 0.;
+                    }
                     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
                 }
             }
             #endif // BASIC_PARAMETERS_H_INCLUDED