rti-theory.c

Go to the documentation of this file.

#include "header.h"

double calc_mse_rti(const gsl_vector *x, void *params)
{
        int i = -1;
        mse_rti_params_struct_type *p = (mse_rti_params_struct_type *) params;
        int nt = p->nt;


        p->alpha = gsl_vector_get(x, 0);
        p->theta = gsl_vector_get(x, 1);
        if (p->alpha <= 0.001) p->alpha = 0.001;
        if (p->theta <= 0.001) p->theta = 0.001;
   
        /* Call the rti_theory function */
        rti_theory(nt, p->spdist, p->samplitude, p->sdelay, p->sduration, p->kappa, p->dfree, 
                   p->alpha, p->theta, p->t, p->p_theory);

        double mse = 0.;

        for (i=1; i<nt; i++) {
                mse += SQR(p->p_model[i] - p->p_theory[i]);
        }
        mse /= nt;

        return mse;
}


void rti_theory(int nt, double spdist, double samplitude, double sdelay, double sduration, double kappa, double dfree, double alpha, double theta, double *t, double *p_theory)
{
        int i;

        double dstar = theta * dfree;
        double ampl = samplitude / (4.0 * PI * alpha * dstar * spdist);

        for (i=0; i<nt; i++) {
                if (t[i] <= sdelay) {
                        p_theory[i] = 0.;
                } else {
                        p_theory[i] = ampl * erfc( spdist /
                                                                        (2.0 * sqrt(dstar * (t[i] - sdelay))) );
                        if (t[i] > sdelay + sduration)
                                p_theory[i] = p_theory[i] - ampl * erfc( spdist /
                                                        (2.0 * sqrt(dstar * (t[i] - (sdelay + sduration)))) );
                }
        }
}