text:
/*
    Licensed to the OpenAirInterface (OAI) Software Alliance under one or more
    contributor license agreements.  See the NOTICE file distributed with
    this work for additional information regarding copyright ownership.
    The OpenAirInterface Software Alliance licenses this file to You under
    the OAI Public License, Version 1.1  (the "License"); you may not use this file
    except in compliance with the License.
    You may obtain a copy of the License at

        http://www.openairinterface.org/?page_id=698

    Unless required by applicable law or agreed to in writing, software
    distributed under the License is distributed on an "AS IS" BASIS,
    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    See the License for the specific language governing permissions and
    limitations under the License.
    -------------------------------------------------------------------------------
    For more information about the OpenAirInterface (OAI) Software Alliance:
        contact@openairinterface.org
*/

#include <math.h>
#include <cblas.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "PHY/TOOLS/tools_defs.h"
#include "sim.h"
#include "scm_corrmat.h"
#include "common/utils/LOG/log.h"
//#define DEBUG_CH

#include "assertions.h"

extern void print_shorts(char *s, __m128i *x);

void fill_channel_desc(channel_desc_t *chan_desc,
                       uint8_t nb_tx,
                       uint8_t nb_rx,
                       uint8_t nb_taps,
                       uint8_t channel_length,
                       double *amps,
                       double *delays,
                       struct complex **R_sqrt,
                       double Td,
                       double sampling_rate,
                       double channel_bandwidth,
                       double ricean_factor,
                       double aoa,
                       double forgetting_factor,
                       double max_Doppler,
                       int32_t channel_offset,
                       double path_loss_dB,
                       uint8_t random_aoa)
{

    uint16_t i, j;
    double delta_tau;

    LOG_I(OCM, "[CHANNEL] Getting new channel descriptor, nb_tx %d, nb_rx %d, nb_taps %d, channel_length %d\n",
          nb_tx, nb_rx, nb_taps, channel_length);

    chan_desc->nb_tx          = nb_tx;
    chan_desc->nb_rx          = nb_rx;
    chan_desc->nb_taps        = nb_taps;
    chan_desc->channel_length = channel_length;
    chan_desc->amps           = amps;
    LOG_D(OCM, "[CHANNEL] Doing delays ...\n");
    if(delays == NULL)
    {
        chan_desc->delays = (double *) malloc(nb_taps * sizeof(double));
        delta_tau = Td / nb_taps;
        for(i = 0; i < nb_taps; i++)
        {
            chan_desc->delays[i] = ((double)i) * delta_tau;
        }
    }
    else
    {
        chan_desc->delays              = delays;
    }
    chan_desc->Td                         = Td;
    chan_desc->sampling_rate              = sampling_rate;
    chan_desc->channel_bandwidth          = channel_bandwidth;
    chan_desc->ricean_factor              = ricean_factor;
    chan_desc->aoa                        = aoa;
    chan_desc->random_aoa                 = random_aoa;
    chan_desc->forgetting_factor          = forgetting_factor;
    chan_desc->channel_offset             = channel_offset;
    chan_desc->path_loss_dB               = path_loss_dB;
    chan_desc->first_run                  = 1;
    chan_desc->ip                         = 0.0;
    chan_desc->max_Doppler                = max_Doppler;
    chan_desc->ch                         = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
    chan_desc->chF                        = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
    chan_desc->a                          = (struct complex **) malloc(nb_taps * sizeof(struct complex *));

    LOG_D(OCM, "[CHANNEL] Filling ch \n");

    for(i = 0; i < nb_tx * nb_rx; i++)
    {
        chan_desc->ch[i] = (struct complex *) malloc(channel_length * sizeof(struct complex));
    }

    for(i = 0; i < nb_tx * nb_rx; i++)
    {
        chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));    // allocate for up to 100 RBs, 12 samples per RB
    }

    LOG_D(OCM, "[CHANNEL] Filling a (nb_taps %d)\n", nb_taps);
    for(i = 0; i < nb_taps; i++)
    {
        LOG_D(OCM, "tap %d (%p,%zu)\n", i, &chan_desc->a[i], nb_tx * nb_rx * sizeof(struct complex));
        chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
    }

    LOG_D(OCM, "[CHANNEL] Doing R_sqrt ...\n");
    if(R_sqrt == NULL)
    {
        chan_desc->R_sqrt         = (struct complex **) calloc(nb_taps, sizeof(struct complex *));
        for(i = 0; i < nb_taps; i++)
        {
            chan_desc->R_sqrt[i]    = (struct complex *) calloc(nb_tx * nb_rx * nb_tx * nb_rx, sizeof(struct complex));
            for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
            {
                chan_desc->R_sqrt[i][j].x = 1.0;
                chan_desc->R_sqrt[i][j].y = 0.0;
            }
        }
    }
    else
    {
        chan_desc->R_sqrt = (struct complex **) calloc(nb_taps, sizeof(struct complex *));
        for(i = 0; i < nb_taps; i++)
        {
            //chan_desc->R_sqrt[i]    = (struct complex*) calloc(nb_tx*nb_rx*nb_tx*nb_rx,sizeof(struct complex));
            //chan_desc->R_sqrt = (struct complex*)&R_sqrt[i][0];
            /* all chan_desc share the same R_sqrt, coming from caller */
            chan_desc->R_sqrt[i] = R_sqrt[0];
        }
    }

    for(i = 0; i < nb_taps; i++)
    {
        for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
        {
            LOG_D(OCM, "Rsqrt[%d][%d] %f %f\n", i, j, chan_desc->R_sqrt[i][j].x, chan_desc->R_sqrt[i][j].y);
        }
    }

    LOG_D(OCM, "[CHANNEL] RF %f\n", chan_desc->ricean_factor);
    for(i = 0; i < chan_desc->nb_taps; i++)
    {
        LOG_D(OCM, "[CHANNEL] tap %d: amp %f, delay %f\n", i, chan_desc->amps[i], chan_desc->delays[i]);
    }

    chan_desc->nb_paths = 10;

    reset_meas(&chan_desc->random_channel);
    reset_meas(&chan_desc->interp_time);
    reset_meas(&chan_desc->interp_freq);
    reset_meas(&chan_desc->convolution);

}

double mbsfn_delays[] = {0, .03, .15, .31, .37, 1.09, 12.490, 12.52, 12.64, 12.80, 12.86, 13.58, 27.49, 27.52, 27.64, 27.80, 27.86, 28.58};
double mbsfn_amps_dB[] = {0, -1.5, -1.4, -3.6, -0.6, -7.0, -10, -11.5, -11.4, -13.6, -10.6, -17.0, -20, -21.5, -21.4, -23.6, -20.6, -27};

double scm_c_delays[] = {0, 0.0125, 0.0250, 0.3625, 0.3750, 0.3875, 0.2500, 0.2625, 0.2750, 1.0375, 1.0500, 1.0625, 2.7250, 2.7375, 2.7500, 4.6000, 4.6125, 4.6250};
double scm_c_amps_dB[] = {0.00, -2.22, -3.98, -1.86, -4.08, -5.84, -1.08, -3.30, -5.06, -9.08, -11.30, -13.06, -15.14, -17.36, -19.12, -20.64, -22.85, -24.62};

double epa_delays[] = { 0, .03, .07, .09, .11, .19, .41};
double epa_amps_dB[] = {0.0, -1.0, -2.0, -3.0, -8.0, -17.2, -20.8};

double eva_delays[] = { 0, .03, .15, .31, .37, .71, 1.09, 1.73, 2.51};
double eva_amps_dB[] = {0.0, -1.5, -1.4, -3.6, -0.6, -9.1, -7.0, -12.0, -16.9};

double etu_delays[] = { 0, .05, .12, .2, .23, .5, 1.6, 2.3, 5.0};
double etu_amps_dB[] = {-1.0, -1.0, -1.0, 0.0, 0.0, 0.0, -3.0, -5.0, -7.0};

double default_amps_lin[] = {0.3868472, 0.3094778, 0.1547389, 0.0773694, 0.0386847, 0.0193424, 0.0096712, 0.0038685};
double default_amp_lin[] = {1};

double ts_shift_delays[] = {0, 1 / 7.68};
double ts_shift_amps[] = {0, 1};

//correlation matrix for a 2x2 channel with full Tx correlation
struct complex R_sqrt_22_corr_tap[16] = {{0.70711, 0}, {0.0, 0.0}, {0.70711, 0}, {0.0, 0.0},
    {0.0, 0.0}, {0.70711, 0}, {0.0, 0.0}, {0.70711, 0},
    {0.70711, 0}, {0.0, 0.0}, {0.70711, 0}, {0.0, 0.0},
    {0.0, 0.0}, {0.70711, 0}, {0.0, 0.0}, {0.70711, 0}
};
struct complex *R_sqrt_22_corr[1]     = {R_sqrt_22_corr_tap};

//correlation matrix for a fully correlated 2x1 channel (h1==h2)
struct complex R_sqrt_21_corr_tap[4]  = {{0.70711, 0}, {0.70711, 0}, {0.70711, 0}, {0.70711, 0}};
struct complex *R_sqrt_21_corr[1]      = {R_sqrt_21_corr_tap};

//correlation matrix for a 2x2 channel with full Tx anti-correlation
struct complex R_sqrt_22_anticorr_tap[16] = {{0.70711, 0}, {0.0, 0.0}, {-0.70711, 0}, {0.0, 0.0},
    {0.0, 0.0}, {0.70711, 0}, {0.0, 0.0}, {-0.70711, 0},
    {-0.70711, 0}, {0.0, 0.0}, {0.70711, 0}, {0.0, 0.0},
    {0.0, 0.0}, {-0.70711, 0}, {0.0, 0.0}, {0.70711, 0}
};
struct complex *R_sqrt_22_anticorr[1]     = {R_sqrt_22_anticorr_tap};

//correlation matrix for a fully anti-correlated 2x1 channel (h1==-h2)
struct complex R_sqrt_21_anticorr_tap[4]  = {{0.70711, 0}, {-0.70711, 0}, {-0.70711, 0}, {0.70711, 0}};
struct complex *R_sqrt_21_anticorr[1]     = {R_sqrt_21_anticorr_tap};

struct complex **R_sqrt_ptr2;

// full correlation matrix in vectorized form for 2x2 channel, where h1 is  perfectly orthogonal to h2

struct complex R_sqrt_22_orthogonal_tap[16] = {{0.70711, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.70711, 0.0}
};
struct complex *R_sqrt_22_orthogonal[1]     = {R_sqrt_22_orthogonal_tap};

// full correlation matrix for TM4 to make orthogonal effective channel

struct complex R_sqrt_22_orth_eff_ch_TM4_prec_real_tap[16] = {{0.70711, 0.0}, {0.0, 0.0}, {0.70711, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.70711, 0.0}, {0.0, 0.0}, {-0.70711, 0.0},
    {0.70711, 0.0}, {0.0, 0.0}, {0.70711, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {-0.70711, 0.0}, {0.0, 0.0}, {0.70711, 0.0}
};
struct complex *R_sqrt_22_orth_eff_ch_TM4_prec_real[1]     = {R_sqrt_22_orth_eff_ch_TM4_prec_real_tap};

struct complex R_sqrt_22_orth_eff_ch_TM4_prec_imag_tap[16] = {{0.70711, 0.0}, {0.0, 0.0}, {0.0, -0.70711}, {0.0, 0.0},
    {0.0, 0.0}, {0.70711, 0.0}, {0.0, 0.0}, {0.0, 0.70711},
    {0.0, -0.70711}, {0.0, 0.0}, {-0.70711, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.0, 0.70711}, {0.0, 0.0}, {-0.70711, 0.0}
};
struct complex *R_sqrt_22_orth_eff_ch_TM4_prec_imag[1]     = {R_sqrt_22_orth_eff_ch_TM4_prec_imag_tap};

//Correlation matrix for EPA channel
struct complex R_sqrt_22_EPA_low_tap[16] = {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}, {0.0, 0.0},
    {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}
};
struct complex *R_sqrt_22_EPA_low[1]     = {R_sqrt_22_EPA_low_tap};

struct complex R_sqrt_22_EPA_high_tap[16] = {{0.7179, 0.0}, {0.4500, 0.0}, {0.4500, 0.0}, {0.2821, 0.0},
    {0.4500, 0.0}, {0.7179, 0.0}, {0.2821, 0.0}, {0.4500, 0.0},
    {0.4500, 0.0}, {0.2821, 0.0}, {0.7179, 0.0}, {0.4500, 0.0},
    {0.2821, 0.0}, {0.4500, 0.0}, {0.4500, 0.0}, {0.7179, 0.0}
};
struct complex *R_sqrt_22_EPA_high[1]     = {R_sqrt_22_EPA_high_tap};

struct complex R_sqrt_22_EPA_medium_tap[16] = {{0.8375, 0.0}, {0.5249, 0.0}, {0.1286, 0.0}, {0.0806, 0.0},
    {0.5249, 0.0}, {0.8375, 0.0}, {0.0806, 0.0}, {0.1286, 0.0},
    {0.1286, 0.0}, {0.0806, 0.0}, {0.8375, 0.0}, {0.5249, 0.0},
    {0.0806, 0.0}, {0.1286, 0.0}, {0.5249, 0.0}, {0.8375, 0.0}
};
struct complex *R_sqrt_22_EPA_medium[1]     = {R_sqrt_22_EPA_medium_tap};

//Rayleigh1_orth_eff_ch_TM4

channel_desc_t *new_channel_desc_scm(uint8_t nb_tx,
                                     uint8_t nb_rx,
                                     SCM_t channel_model,
                                     double sampling_rate,
                                     double channel_bandwidth,
                                     double forgetting_factor,
                                     int32_t channel_offset,
                                     double path_loss_dB)
{

    channel_desc_t *chan_desc = (channel_desc_t *)malloc(sizeof(channel_desc_t));
    uint16_t i, j;
    double sum_amps;
    double aoa, ricean_factor, Td, maxDoppler;
    int channel_length, nb_taps;

    chan_desc->nb_tx                      = nb_tx;
    chan_desc->nb_rx                      = nb_rx;
    chan_desc->sampling_rate              = sampling_rate;
    chan_desc->channel_bandwidth          = channel_bandwidth;
    chan_desc->forgetting_factor          = forgetting_factor;
    chan_desc->channel_offset             = channel_offset;
    chan_desc->path_loss_dB               = path_loss_dB;
    chan_desc->first_run                  = 1;
    chan_desc->ip                                 = 0.0;

    LOG_I(OCM, "Channel Model (inside of new_channel_desc_scm)=%d\n\n", channel_model);

    switch(channel_model)
    {
        case SCM_A:
            LOG_W(OCM, "channel model not yet supported\n");
            free(chan_desc);
            return(NULL);
        case SCM_B:
            LOG_W(OCM, "channel model not yet supported\n");
            free(chan_desc);
            return(NULL);
        case SCM_C:
            chan_desc->nb_taps        = 18;
            chan_desc->Td             = 4.625;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * scm_c_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = scm_c_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }

            chan_desc->R_sqrt  = (struct complex **) malloc(6 * sizeof(struct complex **));
            if(nb_tx == 2 && nb_rx == 2)
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R22_sqrt[i][0];
                }
            }
            else if(nb_tx == 2 && nb_rx == 1)
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R21_sqrt[i][0];
                }
            }
            else if(nb_tx == 1 && nb_rx == 2)
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R12_sqrt[i][0];
                }
            }
            else
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i]    = (struct complex *) malloc(nb_tx * nb_rx * nb_tx * nb_rx * sizeof(struct complex));
                    for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
                    {
                        chan_desc->R_sqrt[i][j].x = 1.0;
                        chan_desc->R_sqrt[i][j].y = 0.0;
                    }
                    LOG_W(OCM, "correlation matrix not implemented for nb_tx==%d and nb_rx==%d, using identity\n", nb_tx, nb_rx);
                }
            }
            break;
        case SCM_D:
            LOG_W(OCM, "This is not the real SCM-D model! It is just SCM-C with an additional Rice factor!\n");
            chan_desc->nb_taps        = 18;
            chan_desc->Td             = 4.625;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * scm_c_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = scm_c_delays;
            chan_desc->ricean_factor  = 0.1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }

            chan_desc->R_sqrt  = (struct complex **) malloc(6 * sizeof(struct complex **));
            if(nb_tx == 2 && nb_rx == 2)
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R22_sqrt[i][0];
                }
            }
            else if(nb_tx == 2 && nb_rx == 1)
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R21_sqrt[i][0];
                }
            }
            else if(nb_tx == 1 && nb_rx == 2)
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R12_sqrt[i][0];
                }
            }
            else
            {
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i]    = (struct complex *) malloc(nb_tx * nb_rx * nb_tx * nb_rx * sizeof(struct complex));
                    for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
                    {
                        chan_desc->R_sqrt[i][j].x = 1.0;
                        chan_desc->R_sqrt[i][j].y = 0.0;
                    }
                    LOG_W(OCM, "correlation matrix not implemented for nb_tx==%d and nb_rx==%d, using identity\n", nb_tx, nb_rx);
                }
            }
            break;
        case EPA:
            chan_desc->nb_taps        = 7;
            chan_desc->Td             = .410;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * epa_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = epa_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }
            if(nb_tx == 2 && nb_rx == 2)
            {
                chan_desc->R_sqrt  = (struct complex **) malloc(6 * sizeof(struct complex **));
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R22_sqrt[i][0];
                }
            }
            else
            {
                chan_desc->R_sqrt         = (struct complex **) malloc(6 * sizeof(struct complex **));
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i]    = (struct complex *) malloc(nb_tx * nb_rx * nb_tx * nb_rx * sizeof(struct complex));
                    for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
                    {
                        chan_desc->R_sqrt[i][j].x = 1.0;
                        chan_desc->R_sqrt[i][j].y = 0.0;
                    }
                    LOG_W(OCM, "correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
                }
            }
            break;
        case EPA_low:
            chan_desc->nb_taps        = 7;
            chan_desc->Td             = .410;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * epa_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = epa_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }
            if(nb_tx == 2 && nb_rx == 2)
            {
                chan_desc->R_sqrt  = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex **));
                for(i = 0; i < chan_desc->nb_taps; i++)
                {
                    chan_desc->R_sqrt[i] = R_sqrt_22_EPA_low[0];
                }
            }
            else
            {
                printf("Correlation matrices are implemented for 2 x 2 only");
            }
            /*  else {
                chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
                for (i = 0; i<6; i++) {
                chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
                for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
                  chan_desc->R_sqrt[i][j].x = 1.0;
                  chan_desc->R_sqrt[i][j].y = 0.0;
                }
                LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
                }
                }*/
            break;
        case EPA_high:
            chan_desc->nb_taps        = 7;
            chan_desc->Td             = .410;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * epa_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = epa_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }
            if(nb_tx == 2 && nb_rx == 2)
            {
                chan_desc->R_sqrt  = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex **));
                for(i = 0; i < chan_desc->nb_taps; i++)
                {
                    chan_desc->R_sqrt[i] = R_sqrt_22_EPA_high[0];
                }
            }
            else
            {
                printf("Correlation matrices are implemented for 2 x 2 only");
            }
            /*  else {
                chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
                for (i = 0; i<6; i++) {
                chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
                for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
                  chan_desc->R_sqrt[i][j].x = 1.0;
                  chan_desc->R_sqrt[i][j].y = 0.0;
                }
                LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
                }
                }*/
            break;
        case EPA_medium:
            chan_desc->nb_taps        = 7;
            chan_desc->Td             = .410;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * epa_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = epa_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }
            if(nb_tx == 2 && nb_rx == 2)
            {
                chan_desc->R_sqrt  = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex **));
                for(i = 0; i < chan_desc->nb_taps; i++)
                {
                    chan_desc->R_sqrt[i] = R_sqrt_22_EPA_medium[0];
                }
            }
            else
            {
                printf("Correlation matrices are implemented for 2 x 2 only");
            }
            /*  else {
                chan_desc->R_sqrt         = (struct complex**) malloc(6*sizeof(struct complex**));
                for (i = 0; i<6; i++) {
                chan_desc->R_sqrt[i]    = (struct complex*) malloc(nb_tx*nb_rx*nb_tx*nb_rx * sizeof(struct complex));
                for (j = 0; j<nb_tx*nb_rx*nb_tx*nb_rx; j+=(nb_tx*nb_rx+1)) {
                  chan_desc->R_sqrt[i][j].x = 1.0;
                  chan_desc->R_sqrt[i][j].y = 0.0;
                }
                LOG_W(OCM,"correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
                }
                }*/
            break;
        case EVA:
            chan_desc->nb_taps        = 9;
            chan_desc->Td             = 2.51;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * eva_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = eva_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }
            if(nb_tx == 2 && nb_rx == 2)
            {
                chan_desc->R_sqrt  = (struct complex **) malloc(6 * sizeof(struct complex **));
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R22_sqrt[i][0];
                }
            }
            else
            {
                chan_desc->R_sqrt         = (struct complex **) malloc(6 * sizeof(struct complex **));
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i]    = (struct complex *) malloc(nb_tx * nb_rx * nb_tx * nb_rx * sizeof(struct complex));
                    for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
                    {
                        chan_desc->R_sqrt[i][j].x = 1.0;
                        chan_desc->R_sqrt[i][j].y = 0.0;
                    }
                    LOG_W(OCM, "correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
                }
            }
            break;
        case ETU:
            chan_desc->nb_taps        = 9;
            chan_desc->Td             = 5.0;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * etu_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = etu_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }
            if(nb_tx == 2 && nb_rx == 2)
            {
                chan_desc->R_sqrt  = (struct complex **) malloc(6 * sizeof(struct complex **));
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i] = (struct complex *) &R22_sqrt[i][0];
                }
            }
            else
            {
                chan_desc->R_sqrt         = (struct complex **) malloc(6 * sizeof(struct complex **));
                for(i = 0; i < 6; i++)
                {
                    chan_desc->R_sqrt[i]    = (struct complex *) malloc(nb_tx * nb_rx * nb_tx * nb_rx * sizeof(struct complex));
                    for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
                    {
                        chan_desc->R_sqrt[i][j].x = 1.0;
                        chan_desc->R_sqrt[i][j].y = 0.0;
                    }
                    LOG_W(OCM, "correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
                }
            }
            break;
        case MBSFN:
            chan_desc->nb_taps        = 18;
            chan_desc->Td             = 28.58;
            chan_desc->channel_length = (int)(2 * chan_desc->sampling_rate * chan_desc->Td + 1 + 2 / (M_PI * M_PI) * log(4 * M_PI * chan_desc->sampling_rate * chan_desc->Td));
            sum_amps = 0;
            chan_desc->amps           = (double *) malloc(chan_desc->nb_taps * sizeof(double));
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i]      = pow(10, .1 * mbsfn_amps_dB[i]);
                sum_amps += chan_desc->amps[i];
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->amps[i] /= sum_amps;
            }
            chan_desc->delays         = mbsfn_delays;
            chan_desc->ricean_factor  = 1;
            chan_desc->aoa            = 0;
            chan_desc->random_aoa     = 0;
            chan_desc->ch             = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->chF            = (struct complex **) malloc(nb_tx * nb_rx * sizeof(struct complex *));
            chan_desc->a              = (struct complex **) malloc(chan_desc->nb_taps * sizeof(struct complex *));
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->ch[i] = (struct complex *) malloc(chan_desc->channel_length * sizeof(struct complex));
            }
            for(i = 0; i < nb_tx * nb_rx; i++)
            {
                chan_desc->chF[i] = (struct complex *) malloc(1200 * sizeof(struct complex));
            }
            for(i = 0; i < chan_desc->nb_taps; i++)
            {
                chan_desc->a[i]         = (struct complex *) malloc(nb_tx * nb_rx * sizeof(struct complex));
            }

            chan_desc->R_sqrt  = (struct complex **) malloc(6 * sizeof(struct complex *));
            for(i = 0; i < 6; i++)
            {
                chan_desc->R_sqrt[i]    = (struct complex *) malloc(nb_tx * nb_rx * nb_tx * nb_rx * sizeof(struct complex));
                for(j = 0; j < nb_tx * nb_rx * nb_tx * nb_rx; j += (nb_tx * nb_rx + 1))
                {
                    chan_desc->R_sqrt[i][j].x = 1.0;
                    chan_desc->R_sqrt[i][j].y = 0.0;
                }
                LOG_W(OCM, "correlation matrix only implemented for nb_tx==2 and nb_rx==2, using identity\n");
            }
            break;
        case Rayleigh8:

            nb_taps = 8;
            Td = 0.8;
            channel_length = (int)11 + 2 * sampling_rate * Td;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 0;

            fill_channel_desc(chan_desc,
                              nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amps_lin,
                              NULL,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rice8:
            nb_taps = 8;
            Td = 0.8;
            channel_length = (int)11 + 2 * sampling_rate * Td;
            ricean_factor = 0.1;
            aoa = 0.7854;
            maxDoppler = 0;

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amps_lin,
                              NULL,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              1);
            break;

        case Rayleigh1:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 0;

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rayleigh1_800:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 800;

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rayleigh1_corr:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 1))
            {
                R_sqrt_ptr2 = R_sqrt_21_corr;
            }
            else if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_corr;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rayleigh1_anticorr:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 1)) //check this
            {
                R_sqrt_ptr2 = R_sqrt_21_anticorr;
            }
            else if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_anticorr;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rice1:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 0.1;
            aoa = 0.7854;
            maxDoppler = 0;

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case AWGN:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 0.0;
            aoa = 0.0;
            maxDoppler = 0;

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            printf("AWGN: ricean_factor %f\n", chan_desc->ricean_factor);

            break;

        case TS_SHIFT:
            nb_taps = 2;
            Td = ts_shift_delays[1];
            channel_length = 10;
            ricean_factor = 0.0;
            aoa = 0.0;
            maxDoppler = 0;

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              ts_shift_amps,
                              ts_shift_delays,
                              NULL,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            printf("TS_SHIFT: ricean_factor %f\n", chan_desc->ricean_factor);

            break;

        case Rice1_corr:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 0.1;
            aoa = .03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 1))
            {
                R_sqrt_ptr2 = R_sqrt_21_corr;
            }
            else if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_corr;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              1);
            break;

        case Rice1_anticorr:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 0.1;
            aoa = .03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 1))
            {
                R_sqrt_ptr2 = R_sqrt_21_anticorr;
            }
            else if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_anticorr;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              1);
            break;

        case Rayleigh1_orthogonal:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = 0.03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_orthogonal;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rayleigh1_orth_eff_ch_TM4_prec_real:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = 0.03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_real;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              1);
            break;

        case Rayleigh1_orth_eff_ch_TM4_prec_imag:
            nb_taps = 1;
            Td = 0;
            channel_length = 1;
            ricean_factor = 1;
            aoa = 0.03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_imag;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc, nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amp_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        case Rayleigh8_orth_eff_ch_TM4_prec_real:

            if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_real;
                //R_sqrt_ptr2 = NULL;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            nb_taps = 8;
            Td = 0.8;
            channel_length = (int)11 + 2 * sampling_rate * Td;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 0;

            fill_channel_desc(chan_desc,
                              nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amps_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);

            break;

        case Rayleigh8_orth_eff_ch_TM4_prec_imag:
            nb_taps = 8;
            Td = 0.8;
            channel_length = (int)11 + 2 * sampling_rate * Td;
            ricean_factor = 1;
            aoa = .03;
            maxDoppler = 0;

            if((nb_tx == 2) && (nb_rx == 2))
            {
                R_sqrt_ptr2 = R_sqrt_22_orth_eff_ch_TM4_prec_imag;
            }
            else
            {
                R_sqrt_ptr2 = NULL;
            }

            fill_channel_desc(chan_desc,
                              nb_tx,
                              nb_rx,
                              nb_taps,
                              channel_length,
                              default_amps_lin,
                              NULL,
                              R_sqrt_ptr2,
                              Td,
                              sampling_rate,
                              channel_bandwidth,
                              ricean_factor,
                              aoa,
                              forgetting_factor,
                              maxDoppler,
                              channel_offset,
                              path_loss_dB,
                              0);
            break;

        default:
            LOG_W(OCM, "channel model not yet supported\n");
            free(chan_desc);
            return(NULL);
    }
    LOG_D(OCM, "[CHANNEL] RF %f\n", chan_desc->ricean_factor);
    for(i = 0; i < chan_desc->nb_taps; i++)
    {
        LOG_D(OCM, "[CHANNEL] tap %d: amp %f, delay %f\n", i, chan_desc->amps[i], chan_desc->delays[i]);
    }

    chan_desc->nb_paths = 10;

    return(chan_desc);
}

void free_channel_desc_scm(channel_desc_t *ch)
{
    // Must be made cleanly, a lot of leaks...
    free(ch);
}

int random_channel(channel_desc_t *desc, uint8_t abstraction_flag)
{

    double s;
    int i, k, l, aarx, aatx;
    struct complex anew[NB_ANTENNAS_TX * NB_ANTENNAS_RX], acorr[NB_ANTENNAS_TX * NB_ANTENNAS_RX];
    struct complex phase, alpha, beta;

    AssertFatal(desc->nb_tx <= NB_ANTENNAS_TX && desc->nb_rx <= NB_ANTENNAS_RX,
                "random_channel.c: Error: temporary buffer for channel not big enough (%d,%d)\n", desc->nb_tx, desc->nb_rx);

    start_meas(&desc->random_channel);
    for(i = 0; i < (int)desc->nb_taps; i++)
    {
        for(aarx = 0; aarx < desc->nb_rx; aarx++)
        {
            for(aatx = 0; aatx < desc->nb_tx; aatx++)
            {

                anew[aarx + (aatx * desc->nb_rx)].x = sqrt(desc->ricean_factor * desc->amps[i] / 2) * gaussdouble(0.0, 1.0);
                anew[aarx + (aatx * desc->nb_rx)].y = sqrt(desc->ricean_factor * desc->amps[i] / 2) * gaussdouble(0.0, 1.0);

                if((i == 0) && (desc->ricean_factor != 1.0))
                {
                    if(desc->random_aoa == 1)
                    {
                        desc->aoa = uniformrandom() * 2 * M_PI;
                    }

                    // this assumes that both RX and TX have linear antenna arrays with lambda/2 antenna spacing.
                    // Furhter it is assumed that the arrays are parallel to each other and that they are far enough apart so
                    // that we can safely assume plane wave propagation.
                    phase.x = cos(M_PI * ((aarx - aatx) * sin(desc->aoa)));
                    phase.y = sin(M_PI * ((aarx - aatx) * sin(desc->aoa)));

                    anew[aarx + (aatx * desc->nb_rx)].x += phase.x * sqrt(1.0 - desc->ricean_factor);
                    anew[aarx + (aatx * desc->nb_rx)].y += phase.y * sqrt(1.0 - desc->ricean_factor);
                }
#ifdef DEBUG_CH
                printf("(%d,%d,%d) %f->(%f,%f) (%f,%f) phase (%f,%f)\n", aarx, aatx, i, desc->amps[i], anew[aarx + (aatx * desc->nb_rx)].x, anew[aarx + (aatx * desc->nb_rx)].y, desc->aoa, desc->ricean_factor, phase.x, phase.y);
#endif
            } //aatx
        } //aarx

        /*
            // for debugging set a=anew;
            for (aarx=0;aarx<desc->nb_rx;aarx++) {
            for (aatx=0;aatx<desc->nb_tx;aatx++) {
            desc->a[i][aarx+(aatx*desc->nb_rx)].x = anew[aarx+(aatx*desc->nb_rx)].x;
            desc->a[i][aarx+(aatx*desc->nb_rx)].y = anew[aarx+(aatx*desc->nb_rx)].y;
            printf("anew(%d,%d) = %f+1j*%f\n",aatx,aarx,anew[aarx+(aatx*desc->nb_rx)].x, anew[aarx+(aatx*desc->nb_rx)].y);
            }
            }
        */
        //apply correlation matrix
        //compute acorr = R_sqrt[i] * anew
        alpha.x = 1.0;
        alpha.y = 0.0;
        beta.x = 0.0;
        beta.y = 0.0;

        cblas_zgemv(CblasRowMajor, CblasNoTrans, desc->nb_tx * desc->nb_rx, desc->nb_tx * desc->nb_rx,
                    (void *) &alpha, (void *) desc->R_sqrt[i / 3], desc->nb_rx * desc->nb_tx,
                    (void *) anew, 1, (void *) &beta, (void *) acorr, 1);

        /*
            for (aarx=0;aarx<desc->nb_rx;aarx++) {
            for (aatx=0;aatx<desc->nb_tx;aatx++) {
            desc->a[i][aarx+(aatx*desc->nb_rx)].x = acorr[aarx+(aatx*desc->nb_rx)].x;
            desc->a[i][aarx+(aatx*desc->nb_rx)].y = acorr[aarx+(aatx*desc->nb_rx)].y;
            printf("tap %d, acorr1(%d,%d) = %f+1j*%f\n",i,aatx,aarx,acorr[aarx+(aatx*desc->nb_rx)].x, acorr[aarx+(aatx*desc->nb_rx)].y);
            }
            }
        */

        if(desc->first_run == 1)
        {
            cblas_zcopy(desc->nb_tx * desc->nb_rx, (void *) acorr, 1, (void *) desc->a[i], 1);
        }
        else
        {
            // a = alpha*acorr+beta*a
            // a = beta*a
            // a = a+alpha*acorr
            alpha.x = sqrt(1 - desc->forgetting_factor);
            alpha.y = 0;
            beta.x = sqrt(desc->forgetting_factor);
            beta.y = 0;
            cblas_zscal(desc->nb_tx * desc->nb_rx, (void *) &beta, (void *) desc->a[i], 1);
            cblas_zaxpy(desc->nb_tx * desc->nb_rx, (void *) &alpha, (void *) acorr, 1, (void *) desc->a[i], 1);

            //  desc->a[i][aarx+(aatx*desc->nb_rx)].x = (sqrt(desc->forgetting_factor)*desc->a[i][aarx+(aatx*desc->nb_rx)].x) + sqrt(1-desc->forgetting_factor)*anew.x;
            //  desc->a[i][aarx+(aatx*desc->nb_rx)].y = (sqrt(desc->forgetting_factor)*desc->a[i][aarx+(aatx*desc->nb_rx)].y) + sqrt(1-desc->forgetting_factor)*anew.y;
        }

        /*
            for (aarx=0;aarx<desc->nb_rx;aarx++) {
            for (aatx=0;aatx<desc->nb_tx;aatx++) {
            //desc->a[i][aarx+(aatx*desc->nb_rx)].x = acorr[aarx+(aatx*desc->nb_rx)].x;
            //desc->a[i][aarx+(aatx*desc->nb_rx)].y = acorr[aarx+(aatx*desc->nb_rx)].y;
            printf("tap %d, a(%d,%d) = %f+1j*%f\n",i,aatx,aarx,desc->a[i][aarx+(aatx*desc->nb_rx)].x, desc->a[i][aarx+(aatx*desc->nb_rx)].y);
            }
            }
        */

    } //nb_taps
    stop_meas(&desc->random_channel);

    //memset((void *)desc->ch[aarx+(aatx*desc->nb_rx)],0,(int)(desc->channel_length)*sizeof(struct complex));

    if(abstraction_flag == 0)
    {
        start_meas(&desc->interp_time);

        for(aarx = 0; aarx < desc->nb_rx; aarx++)
        {
            for(aatx = 0; aatx < desc->nb_tx; aatx++)
            {
                if(desc->channel_length == 1)
                {
                    desc->ch[aarx + (aatx * desc->nb_rx)][0].x = desc->a[0][aarx + (aatx * desc->nb_rx)].x;
                    desc->ch[aarx + (aatx * desc->nb_rx)][0].y = desc->a[0][aarx + (aatx * desc->nb_rx)].y;
                }
                else
                {

                    for(k = 0; k < (int)desc->channel_length; k++)
                    {
                        desc->ch[aarx + (aatx * desc->nb_rx)][k].x = 0.0;
                        desc->ch[aarx + (aatx * desc->nb_rx)][k].y = 0.0;

                        for(l = 0; l < desc->nb_taps; l++)
                        {
                            if((k - (desc->delays[l]*desc->sampling_rate) - NB_SAMPLES_CHANNEL_OFFSET) == 0)
                            {
                                s = 1.0;
                            }
                            else
                                s = sin(M_PI * (k - (desc->delays[l] * desc->sampling_rate) - NB_SAMPLES_CHANNEL_OFFSET)) /
                                    (M_PI * (k - (desc->delays[l] * desc->sampling_rate) - NB_SAMPLES_CHANNEL_OFFSET));

                            desc->ch[aarx + (aatx * desc->nb_rx)][k].x += s * desc->a[l][aarx + (aatx * desc->nb_rx)].x;
                            desc->ch[aarx + (aatx * desc->nb_rx)][k].y += s * desc->a[l][aarx + (aatx * desc->nb_rx)].y;
                            //    printf("l %d : desc->ch.x %f\n",l,desc->a[l][aarx+(aatx*desc->nb_rx)].x);

                        } //nb_taps

#ifdef DEBUG_CH
                        k = 0;
                        printf("(%d,%d,%d)->(%f,%f)\n", k, aarx, aatx, desc->ch[aarx + (aatx * desc->nb_rx)][k].x, desc->ch[aarx + (aatx * desc->nb_rx)][k].y);
#endif
                    }
                } //channel_length
            } //aatx
        } //aarx
        stop_meas(&desc->interp_time);
    }

    if(desc->first_run == 1)
    {
        desc->first_run = 0;
    }

    return (0);
}

double N_RB2sampling_rate(uint16_t N_RB)
{
    double sampling_rate;
    switch(N_RB)
    {
        case 6:
            sampling_rate = 1.92;
            break;

        case 25:
            sampling_rate = 7.68;
            break;

        case 50:
            sampling_rate = 15.36;
            break;

        case 100:
            sampling_rate = 30.72;
            break;

        default:
            AssertFatal(1 == 0, "Unknown N_PRB %d", N_RB);
    }

    return(sampling_rate);
}

double N_RB2channel_bandwidth(uint16_t N_RB)
{
    double channel_bandwidth;
    switch(N_RB)
    {
        case 6:
            channel_bandwidth = 1.25;
            break;

        case 25:
            channel_bandwidth = 5.00;
            break;

        case 50:
            channel_bandwidth = 10.00;
            break;

        case 100:
            channel_bandwidth = 20.00;
            break;

        default:
            LOG_E(PHY, "Unknown N_PRB\n");
            return(-1);
    }
    return(channel_bandwidth);
}

#ifdef RANDOM_CHANNEL_MAIN
#define sampling_rate 5.0
#define Td 2.0
main(int argc, char **argv)
{

    double amps[8] = {.8, .2, .1, .04, .02, .01, .005};
    struct complex ch[(int)(1 + 2 * sampling_rate * Td)], phase;
    int i;

    randominit();
    phase.x = 1.0;
    phase.y = 0;
    random_channel(amps, Td, 8, sampling_rate, ch, (double)1.0, &phase);
    /*
        for (i=0;i<(11+2*sampling_rate*Td);i++){
        printf("%f + sqrt(-1)*%f\n",ch[i].x,ch[i].y);
        }
    */
}

#endif


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <opencv2\core.hpp>
#include <gsl\gsl>

namespace mage
{
    namespace clouds
    {
        void ComputeCharacteristics(
            gsl::span<const cv::Point3f> points,
            std::vector<cv::Vec3f>& normals,
            gsl::span<const size_t> mapIds,
            std::vector<float>& effectiveDissimilarities,
            std::vector<float>& distanceScores,
            std::vector<float>& homogeneityScores,
            const size_t maxN);

        void MollifyNormals(
            gsl::span<const cv::Point3f> points,
            gsl::span<cv::Vec3f> normals,
            const float normalScale,
            const float spatialScale,
            const size_t iterations,
            const size_t maxN);

        void RepositionPointSets(
            gsl::span<cv::Point3f> toReposition,
            gsl::span<cv::Vec3f> normals,
            const float stepi,
            const size_t iterations,
            const float lambda,
            const float solverSigmaS,
            const float mollSigmaS,
            const float mollSigmaN,
            const size_t maxN);

        struct Knn final
        {
            Knn(gsl::span<const cv::Point3f> points);
            Knn(gsl::span<const cv::Vec3f> points);
            ~Knn();

            Knn& operator=(Knn&& other);

            void Query(const cv::Vec3f& point, const size_t N, std::vector<size_t>& indexes) const;

        private:
            struct Impl;
            std::unique_ptr<Impl> m_impl;
        };
    }
}

ext: h

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#pragma once

#include "Availability.h"
#include "Library.h"
#include "OctreeTileID.h"
#include "TileAvailabilityFlags.h"

#include <gsl/span>

#include <cstddef>
#include <memory>
#include <vector>

namespace CesiumGeometry {

class CESIUMGEOMETRY_API OctreeAvailability final {
public:
  /**
   * @brief Constructs a new instance.
   *
   * @param subtreeLevels The number of levels in each subtree.
   * @param maximumLevel The index of the maximum level in this tileset.
   */
  OctreeAvailability(uint32_t subtreeLevels, uint32_t maximumLevel) noexcept;

  /**
   * @brief Determines the currently known availability status of the given
   * tile.
   *
   * @param tileID The {@link CesiumGeometry::OctreeTileID} for the tile.
   *
   * @return The {@link TileAvailabilityFlags} for this tile encoded into a
   * uint8_t.
   */
  uint8_t computeAvailability(const OctreeTileID& tileID) const noexcept;

  /**
   * @brief Attempts to add an availability subtree into the existing overall
   * availability tree.
   *
   * @param tileID The {@link CesiumGeometry::OctreeTileID} for the tile.
   * @param newSubtree The {@link CesiumGeometry::AvailabilitySubtree} to add.
   *
   * @return Whether the insertion was successful.
   */
  bool addSubtree(
      const OctreeTileID& tileID,
      AvailabilitySubtree&& newSubtree) noexcept;

  /**
   * @brief Determines the currently known availability status of the given
   * tile.
   *
   * Priming with a known parent subtree node avoids the need to traverse the
   * entire availability tree so far. The node must have a loaded subtree
   *
   * @param tileID The tile ID to get the availability for.
   * @param pNode The subtree node to look for the tileID in. The tileID should
   * be within this subtree node.
   *
   * @return The {@link TileAvailabilityFlags} for this tile encoded into a
   * uint8_t.
   */
  uint8_t computeAvailability(
      const OctreeTileID& tileID,
      const AvailabilityNode* pNode) const noexcept;

  /**
   * @brief Attempts to add a child subtree node onto the given parent node.
   *
   * Priming with a known parent subtree node avoids the need to traverse the
   * entire availability tree so far. If the parent node is nullptr and the
   * tile ID indicates this is the root tile, the subtree will be attached to
   * the root.
   *
   * @param tileID The root tile's ID of the subtree we are trying to add.
   * @param pParentNode The parent subtree node. The tileID should fall exactly
   * at the end of this parent subtree.
   *
   * @return The newly created node if the insertion was successful, nullptr
   * otherwise.
   */
  AvailabilityNode*
  addNode(const OctreeTileID& tileID, AvailabilityNode* pParentNode) noexcept;

  /**
   * @brief Attempts to add a loaded subtree onto the given node.
   *
   * The node must have been created earlier from a call to addNode.
   *
   * @param pNode The node on which to add the subtree.
   * @param newSubtree The new subtree to add.
   *
   * @return Whether the insertion was successful.
   */
  bool addLoadedSubtree(
      AvailabilityNode* pNode,
      AvailabilitySubtree&& newSubtree) noexcept;
  /**
   * @brief Find the child node index corresponding to this tile ID and parent
   * node.
   *
   * Attempts to find the child node for the tile with the given ID and parent
   * node. The parent node is used to speed up the search significantly. Note
   * that if the given tile ID does not correspond exactly to an immediate
   * child node of the parent node, nullptr will be returned. If a tileID
   * outside the given parent node's subtree is given, an incorrect child index
   * may be returned.
   *
   * @param tileID The tile ID of the child node we are looking for.
   * @param pParentNode The immediate parent to the child node we are looking
   * for.
   * @return The child node index if found, std::nullopt otherwise.
   */
  std::optional<uint32_t> findChildNodeIndex(
      const OctreeTileID& tileID,
      const AvailabilityNode* pParentNode) const;

  /**
   * @brief Find the child node corresponding to this tile ID and parent node.
   *
   * Attempts to find the child node for the tile with the given ID and parent
   * node. The parent node is used to speed up the search significantly. Note
   * that if the given tile ID does not correspond exactly to an immediate
   * child node of the parent node, nullptr will be returned. If a tileID
   * outside the given parent node's subtree is given, an incorrect child index
   * may be returned.
   *
   * @param tileID The tile ID of the child node we are looking for.
   * @param pParentNode The immediate parent to the child node we are looking
   * for.
   * @return The child node if found, nullptr otherwise.
   */
  AvailabilityNode* findChildNode(
      const OctreeTileID& tileID,
      AvailabilityNode* pParentNode) const;

  /**
   * @brief Gets the number of levels in each subtree.
   */
  constexpr inline uint32_t getSubtreeLevels() const noexcept {
    return this->_subtreeLevels;
  }

  /**
   * @brief Gets the index of the maximum level in this implicit tileset.
   */
  constexpr inline uint32_t getMaximumLevel() const noexcept {
    return this->_maximumLevel;
  }

  /**
   * @brief Gets a pointer to the root subtree node of this implicit tileset.
   */
  AvailabilityNode* getRootNode() noexcept { return this->_pRoot.get(); }

private:
  uint32_t _subtreeLevels;
  uint32_t _maximumLevel;
  uint32_t _maximumChildrenSubtrees;
  std::unique_ptr<AvailabilityNode> _pRoot;
};

} // namespace CesiumGeometry

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
static const char help[] =
"Solves obstacle problem in 2D using SNESVI.  Option prefix -obs_.\n"
"The obstacle problem is a free boundary problem for the Poisson equation\n"
"in which the solution u(x,y) is constrained to be above the obstacle psi(x,y):\n"
"    - Lap u = f,  u >= psi.\n"
"Equivalently it is a variational inequality (VI), complementarity problem\n"
"(CP), or an inequality-constrained minimization.  The example here is\n"
"on the square (-2,2)^2 and has known exact solution.  Because of the\n"
"constraint, the problem is nonlinear but the code reuses the residual and\n"
"Jacobian evaluation code for the Poisson equation in ch6/.\n\n";

#include <petsc.h>
#include "../ch6/poissonfunctions.h"

// z = psi(x,y) is the hemispherical obstacle, but made C^1 with "skirt" at r=r0
PetscReal psi(PetscReal x, PetscReal y) {
    const PetscReal  r = x * x + y * y,
                     r0 = 0.9,
                     psi0 = PetscSqrtReal(1.0 - r0*r0),
                     dpsi0 = - r0 / psi0;
    if (r <= r0) {
        return PetscSqrtReal(1.0 - r);
    } else {
        return psi0 + dpsi0 * (r - r0);
    }
}

/*  This exact solution solves a 1D radial free-boundary problem for the
Laplace equation, on the interval 0 < r < 2, with hemispherical obstacle
    psi(r) =  / sqrt(1 - r^2),  r < 1
              \ -1,             otherwise
The Laplace equation applies where u(r) > psi(r),
    u''(r) + r^-1 u'(r) = 0
with boundary conditions including free b.c.s at an unknown location r = a:
    u(a) = psi(a),  u'(a) = psi'(a),  u(2) = 0
The solution is  u(r) = - A log(r) + B   on  r > a.  The boundary conditions
can then be reduced to a root-finding problem for a:
    a^2 (log(2) - log(a)) = 1 - a^2
The solution is a = 0.697965148223374 (giving residual 1.5e-15).  Then
A = a^2*(1-a^2)^(-0.5) and B = A*log(2) are as given below in the code.  */
PetscReal u_exact(PetscReal x, PetscReal y) {
    const PetscReal afree = 0.697965148223374,
                    A     = 0.680259411891719,
                    B     = 0.471519893402112;
    PetscReal       r;
    r = PetscSqrtReal(x * x + y * y);
    return (r <= afree) ? psi(x,y)  // active set; on the obstacle
                        : - A * PetscLogReal(r) + B; // solves laplace eqn
}

// boundary conditions from exact solution
PetscReal g_fcn(PetscReal x, PetscReal y, PetscReal z, void *ctx) {
    return u_exact(x,y);
}

// we solve Laplace's equation with f = 0
PetscReal zero(PetscReal x, PetscReal y, PetscReal z, void *ctx) {
    return 0.0;
}

extern PetscErrorCode FormUExact(DMDALocalInfo*, Vec);
extern PetscErrorCode GetActiveSet(SNES, DMDALocalInfo*, Vec, Vec,
                                   PetscInt*, PetscReal*);
extern PetscErrorCode FormBounds(SNES, Vec, Vec);

int main(int argc,char **argv) {
  PetscErrorCode ierr;
  DM                  da, da_after;
  SNES                snes;
  KSP                 ksp;
  Vec                 u_initial, u, u_exact, Xl, Xu;
  PoissonCtx          user;
  const PetscReal     aexact = 0.697965148223374;
  SNESConvergedReason reason;
  PetscInt            snesit, kspit;
  PetscReal           error1,errorinf,actarea,exactarea,areaerr;
  DMDALocalInfo       info;
  char                dumpname[256] = "dump.dat";
  PetscBool           dumpbinary = PETSC_FALSE;

  ierr = PetscInitialize(&argc,&argv,NULL,help); if (ierr) return ierr;

  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"obs_","options to obstacle","");CHKERRQ(ierr);
  ierr = PetscOptionsString("-dump_binary",
           "filename for saving solution AND OBSTACLE in PETSc binary format",
           "obstacle.c",dumpname,dumpname,sizeof(dumpname),&dumpbinary); CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  ierr = DMDACreate2d(PETSC_COMM_WORLD,
      DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
      3,3,                       // override with -da_grid_x,_y
      PETSC_DECIDE,PETSC_DECIDE, // num of procs in each dim
      1,1,NULL,NULL,             // dof = 1 and stencil width = 1
      &da);CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(da,-2.0,2.0,-2.0,2.0,-1.0,-1.0);CHKERRQ(ierr);

  user.cx = 1.0;
  user.cy = 1.0;
  user.cz = 1.0;
  user.g_bdry = &g_fcn;
  user.f_rhs = &zero;
  user.addctx = NULL;
  ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da);CHKERRQ(ierr);
  ierr = SNESSetApplicationContext(snes,&user);CHKERRQ(ierr);

  // set the SNES type to a variational inequality (VI) solver of reduced-space
  // (RS) type
  ierr = SNESSetType(snes,SNESVINEWTONRSLS);CHKERRQ(ierr);
  ierr = SNESVISetComputeVariableBounds(snes,&FormBounds);CHKERRQ(ierr);

  // reuse residual and jacobian from ch6/
  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)Poisson2DFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDASNESSetJacobianLocal(da,
             (DMDASNESJacobian)Poisson2DJacobianLocal,&user); CHKERRQ(ierr);
  ierr = SNESGetKSP(snes,&ksp); CHKERRQ(ierr);
  ierr = KSPSetType(ksp,KSPCG); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);

  // initial iterate is zero for simplicity
  ierr = DMCreateGlobalVector(da,&u_initial);CHKERRQ(ierr);
  ierr = VecSet(u_initial,0.0); CHKERRQ(ierr);

  /* solve and get solution, DM after solution*/
  ierr = SNESSolve(snes,NULL,u_initial);CHKERRQ(ierr);
  ierr = VecDestroy(&u_initial); CHKERRQ(ierr);
  ierr = DMDestroy(&da); CHKERRQ(ierr);
  ierr = SNESGetDM(snes,&da_after); CHKERRQ(ierr);
  ierr = SNESGetSolution(snes,&u); CHKERRQ(ierr); /* do not destroy u */
  ierr = DMDAGetLocalInfo(da_after,&info); CHKERRQ(ierr);
  ierr = VecDuplicate(u,&Xl); CHKERRQ(ierr);
  ierr = VecDuplicate(u,&Xu); CHKERRQ(ierr);
  ierr = FormBounds(snes,Xl,Xu); CHKERRQ(ierr);

  /* save solution to binary file if requested */
  if (dumpbinary) {
      PetscViewer dumpviewer;
      ierr = PetscPrintf(PETSC_COMM_WORLD,
               "writing u,psi in binary format to %s ...\n",dumpname); CHKERRQ(ierr);
      ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,dumpname,FILE_MODE_WRITE,&dumpviewer); CHKERRQ(ierr);
      ierr = VecView(u,dumpviewer); CHKERRQ(ierr);
      ierr = VecView(Xl,dumpviewer); CHKERRQ(ierr);
      ierr = PetscViewerDestroy(&dumpviewer); CHKERRQ(ierr);
  }

  /* compute final performance measures */
  ierr = SNESGetConvergedReason(snes,&reason); CHKERRQ(ierr);
  if (reason <= 0) {
      ierr = PetscPrintf(PETSC_COMM_WORLD,
          "WARNING: SNES not converged ... use -snes_converged_reason to check\n"); CHKERRQ(ierr);
  }
  ierr = SNESGetIterationNumber(snes,&snesit); CHKERRQ(ierr);
  ierr = SNESGetKSP(snes,&ksp); CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&kspit); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "done on %d x %d grid ... %s, SNES iters = %d, last KSP iters = %d\n",
      info.mx,info.my,SNESConvergedReasons[reason],snesit,kspit); CHKERRQ(ierr);

  /* compare to exact */
  ierr = GetActiveSet(snes,&info,u,Xl,NULL,&actarea); CHKERRQ(ierr);
  exactarea = PETSC_PI * aexact * aexact;
  areaerr = PetscAbsReal(actarea - exactarea) / exactarea;
  ierr = VecDuplicate(u,&u_exact); CHKERRQ(ierr);
  ierr = FormUExact(&info,u_exact); CHKERRQ(ierr);
  ierr = VecAXPY(u,-1.0,u_exact); CHKERRQ(ierr); /* u <- u - u_exact */
  ierr = VecNorm(u,NORM_1,&error1); CHKERRQ(ierr);
  error1 /= (PetscReal)info.mx * (PetscReal)info.my;
  ierr = VecNorm(u,NORM_INFINITY,&errorinf); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,
      "errors: av |u-uexact| = %.3e, |u-uexact|_inf = %.3e, active area error = %.3f%%\n",
      error1,errorinf,100.0*areaerr); CHKERRQ(ierr);

  VecDestroy(&u_exact);  VecDestroy(&Xl);  VecDestroy(&Xu);
  SNESDestroy(&snes);
  return PetscFinalize();
}

PetscErrorCode GetActiveSet(SNES snes, DMDALocalInfo *info, Vec u, Vec Xl,
                            PetscInt *act, PetscReal *actarea) {
  PetscErrorCode ierr;
  Vec              F;
  const PetscReal  *au, *aXl, *aF, zerotol = 1.0e-8;  // see petsc/src/snes/impls/vi/vi.c for value
  PetscReal        dx, dy;
  PetscInt         i,n,lact,gact;

  dx = 4.0 / (PetscReal)(info->mx-1);
  dy = 4.0 / (PetscReal)(info->my-1);
  ierr = VecGetLocalSize(u,&n);CHKERRQ(ierr);
  ierr = VecGetArrayRead(u,&au); CHKERRQ(ierr);
  ierr = VecGetArrayRead(Xl,&aXl); CHKERRQ(ierr);
  ierr = SNESGetFunction(snes,&F,NULL,NULL); CHKERRQ(ierr); /* do not destroy F */
  ierr = VecGetArrayRead(F,&aF); CHKERRQ(ierr);
  lact = 0;
  for (i=0; i<n; i++) {
    if ((au[i] <= aXl[i] + zerotol) && (aF[i] > 0.0))
        lact++;
  }
  ierr = VecRestoreArrayRead(u,&au);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(Xl,&aXl);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(F,&aF); CHKERRQ(ierr);
  ierr = MPI_Allreduce(&lact,&gact,1,MPIU_INT,MPIU_SUM,PetscObjectComm((PetscObject)snes));CHKERRQ(ierr);
  if (act) {
      *act = gact;
  }
  if (actarea) {
      *actarea = dx * dy * gact;
  }
  return 0;
}

PetscErrorCode FormUExact(DMDALocalInfo *info, Vec u) {
  PetscErrorCode ierr;
  PetscInt   i,j;
  PetscReal  **au, dx, dy, x, y;
  dx = 4.0 / (PetscReal)(info->mx-1);
  dy = 4.0 / (PetscReal)(info->my-1);
  ierr = DMDAVecGetArray(info->da, u, &au);CHKERRQ(ierr);
  for (j=info->ys; j<info->ys+info->ym; j++) {
    y = -2.0 + j * dy;
    for (i=info->xs; i<info->xs+info->xm; i++) {
      x = -2.0 + i * dx;
      au[j][i] = u_exact(x,y);
    }
  }
  ierr = DMDAVecRestoreArray(info->da, u, &au);CHKERRQ(ierr);
  return 0;
}

//STARTBOUNDS
// for call-back: tell SNESVI we want  psi <= u < +infinity
PetscErrorCode FormBounds(SNES snes, Vec Xl, Vec Xu) {
  PetscErrorCode ierr;
  DM             da;
  DMDALocalInfo  info;
  PetscInt       i, j;
  PetscReal      **aXl, dx, dy, x, y;
  ierr = SNESGetDM(snes,&da);CHKERRQ(ierr);
  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  dx = 4.0 / (PetscReal)(info.mx-1);
  dy = 4.0 / (PetscReal)(info.my-1);
  ierr = DMDAVecGetArray(da, Xl, &aXl);CHKERRQ(ierr);
  for (j=info.ys; j<info.ys+info.ym; j++) {
    y = -2.0 + j * dy;
    for (i=info.xs; i<info.xs+info.xm; i++) {
      x = -2.0 + i * dx;
      aXl[j][i] = psi(x,y);
    }
  }
  ierr = DMDAVecRestoreArray(da, Xl, &aXl);CHKERRQ(ierr);
  ierr = VecSet(Xu,PETSC_INFINITY);CHKERRQ(ierr);
  return 0;
}
//ENDBOUNDS

ext: c


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>

void test_matmul(){
    int i, j;

    gsl_matrix * m1 = gsl_matrix_alloc(10, 4);
    gsl_matrix * m2 = gsl_matrix_alloc(10, 4);
    for(i=0; i<10; i++){
        for(j=0; j<4; j++){
            gsl_matrix_set(m1, i, j, rand() % 10);
        }
    }

    for(i=0; i<10; i++){
        for(j=0; j<4; j++){
            gsl_matrix_set(m2, i, j, rand() % 10);
        }
    }

    printf("m1 before matmul:\n");
    for (i = 0; i < 10; i++){
        for (j = 0; j < 4; j++){
            printf ("%g  ", gsl_matrix_get (m1, i, j));
        }
        printf("\n");
    }

    gsl_matrix * result = gsl_matrix_alloc(4, 4);

    gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, m1, m2, 0.0, result);

    printf("m1 after matmul:\n");
    for (i = 0; i < 10; i++){
        for (j = 0; j < 4; j++){
            printf ("%g  ", gsl_matrix_get (m1, i, j));
        }
        printf("\n");
    }

    printf("matmul result:\n");
    for (i = 0; i < 4; i++){
        for (j = 0; j < 4; j++){
            printf ("%g ", gsl_matrix_get (result, i, j));
        }
        printf("\n");
    }

    // gsl_matrix_fprintf(stdout, result, "%g");

    gsl_matrix_free(m1);
    gsl_matrix_free(m2);
    gsl_matrix_free(result);

}

void test_mat_inverse(){
    int i, j;
    double m1_data[] = {4, 7,
                        2, 6};
    gsl_matrix_view m1 = gsl_matrix_view_array(m1_data, 2, 2);
    int s;
    gsl_permutation * p1 = gsl_permutation_alloc(2);
    gsl_matrix * m1_inv = gsl_matrix_alloc(2, 2);
    gsl_linalg_LU_decomp(&m1.matrix, p1, &s);
    gsl_linalg_LU_invert(&m1.matrix, p1, m1_inv);
    printf("Inverse of 2x2 matrix:\n");
    for (i = 0; i < 2; i++){
        for (j = 0; j < 2; j++){
            printf ("%g ", gsl_matrix_get (m1_inv, i, j));
        }
        printf("\n");
    }

    double m2_data[] = {7, 2, 1,
                        0, 3, -1,
                        -3, 4, -2};
    gsl_matrix_view m2 = gsl_matrix_view_array(m2_data, 3, 3);
    gsl_permutation * p2 = gsl_permutation_alloc(3);
    gsl_matrix * m2_inv = gsl_matrix_alloc(3, 3);
    gsl_linalg_LU_decomp(&m2.matrix, p2, &s);
    gsl_linalg_LU_invert(&m2.matrix, p2, m2_inv);
    printf("Inverse of 3x3 matrix:\n");
    for (i = 0; i < 3; i++){
        for (j = 0; j < 3; j++){
            printf ("%g ", gsl_matrix_get (m2_inv, i, j));
        }
        printf("\n");
    }

    double m3_data[] = {6, 4, 2, 8,
                        1, 3, 5, 7,
                        2, 3, 5, 6,
                        4, 2, 6, 3};
    gsl_matrix_view m3 = gsl_matrix_view_array(m3_data, 4, 4);
    gsl_permutation * p3 = gsl_permutation_alloc(4);
    gsl_matrix * m3_inv = gsl_matrix_alloc(4, 4);
    gsl_linalg_LU_decomp(&m3.matrix, p3, &s);
    gsl_linalg_LU_invert(&m3.matrix, p3, m3_inv);
    printf("Inverse of 4x4 matrix:\n");
    for (i = 0; i < 4; i++){
        for (j = 0; j < 4; j++){
            printf ("%g ", gsl_matrix_get (m3_inv, i, j));
        }
        printf("\n");
    }

    // gsl_matrix_fprintf(stdout, m3_inv, "%g");

}

int main(){
    test_matmul();
    //test_mat_inverse();
}


text:
/* rng/fishman2x.c
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
 * This generator is taken from
 *
 * Donald E. Knuth
 * The Art of Computer Programming
 * Volume 2
 * Third Edition
 * Addison-Wesley
 * Page 108
 *
 * It is called "Fishman - L'Ecuyer"
 *
 * This implementation copyright (C) 2001 Carlo Perassi
 * and (C) 2003 Heiko Bauke.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* Fishman */
#define AAA_F 48271UL
#define MMM_F 0x7fffffffUL      /* 2 ^ 31 - 1 */
#define QQQ_F 44488UL
#define RRR_F 3399UL

/* L'Ecuyer */
#define AAA_L 40692UL
#define MMM_L 0x7fffff07UL      /* 2 ^ 31 - 249 */
#define QQQ_L 52774UL
#define RRR_L 3791UL

static inline unsigned long int ran_get (void *vstate);
static double ran_get_double (void *vstate);
static void ran_set (void *state, unsigned long int s);

typedef struct
{
  unsigned long int x;
  unsigned long int y;
  unsigned long int z;
}
ran_state_t;

static inline unsigned long int
ran_get (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  long int y, r;

  r = RRR_F * (state->x / QQQ_F);
  y = AAA_F * (state->x % QQQ_F) - r;
  if (y < 0)
    y += MMM_F;
  state->x = y;

  r = RRR_L * (state->y / QQQ_L);
  y = AAA_L * (state->y % QQQ_L) - r;
  if (y < 0)
    y += MMM_L;
  state->y = y;

  state->z = (state->x > state->y) ? (state->x - state->y) :
    MMM_F + state->x - state->y;

  return state->z;
}

static double
ran_get_double (void *vstate)
{
  ran_state_t *state = (ran_state_t *) vstate;

  return ran_get (state) / 2147483647.0;
}

static void
ran_set (void *vstate, unsigned long int s)
{
  ran_state_t *state = (ran_state_t *) vstate;

  if ((s % MMM_F) == 0 || (s % MMM_L) == 0)
    s = 1;                      /* default seed is 1 */

  state->x = s % MMM_F;
  state->y = s % MMM_L;
  state->z = (state->x > state->y) ? (state->x - state->y) :
    MMM_F + state->x - state->y;

  return;
}

static const gsl_rng_type ran_type = {
  "fishman2x",                  /* name */
  MMM_F - 1,                    /* RAND_MAX */
  0,                            /* RAND_MIN */
  sizeof (ran_state_t),
  &ran_set,
  &ran_get,
  &ran_get_double
};

const gsl_rng_type *gsl_rng_fishman2x = &ran_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpfr.h>
#include <petsc.h>

#include "../ellipsoid/ellipsoid.h"
#include "../ellipsoid/ellSolv.h"
#include "../constants.h"
#include "testFunctions.h"

#undef __FUNCT__
#define __FUNCT__ "main"
PetscErrorCode main(int argc, char **argv)
{

  PetscErrorCode ierr;
  EllipsoidalSystem e;
  PetscInt prec = 64;
  PetscReal a = 3.0;
  PetscReal b = 2.0;
  PetscReal c = 1.0;
  double doubleEnp;
  mpfr_t x;
  mpfr_t mpfrEnp;
  PetscFunctionBeginUser;

  ierr = PetscInitialize(&argc, &argv, NULL, NULL);CHKERRQ(ierr);
  ierr = PetscLogDefaultBegin();CHKERRQ(ierr);

  mpfr_inits(mpfrEnp, x, NULL);

  ierr = initEllipsoidalSystem(&e, a, b, c, prec);CHKERRQ(ierr);

  mpfr_set_d(x, 1.0, MPFR_RNDN);
  //calculate with double
  ierr = calcLame(&e, 3, 4, 1.0, 1, 1, &doubleEnp);
  //calculate with high precision
  ierr = CalcLameMPFR(&e, 3, 4, x, 1, 1, &mpfrEnp);

  printf("Enp(double): %18.18f\n", doubleEnp);
  printf("Enp(mpfr):   %18.18f\n", mpfr_get_d(mpfrEnp, MPFR_RNDN));

  mpfr_clears(mpfrEnp, x, NULL);
  ierr = PetscFinalize();CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

max_stars_repo_name: tom-klotz/ellipsoid-solvation

max_issues_repo_name: tom-klotz/ellipsoid-solvation

max_forks_repo_name: tom-klotz/ellipsoid-solvation