E_1DS.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include "../ascot5.h"
#include "../error.h"
#include "../print.h"
#include "../B_field.h"
#include "../spline/interp.h"
#include "E_1DS.h"
 
int E_1DS_init_offload(E_1DS_offload_data* offload_data, real** offload_array) {
 
    /* Spline initialization. */
    int err = 0;
    real* dV_drho = (real*) malloc(NSIZE_COMP1D * offload_data->n_rho
                                   * sizeof(real));
    err += interp1Dcomp_init_coeff(dV_drho, *offload_array,
                                   offload_data->n_rho,
                                   NATURALBC,
                                   offload_data->rho_min,
                                   offload_data->rho_max);
    if(err) {
        print_err("Error: Failed to initialize splines.\n");
        return err;
    }
 
    /* Free the offload array and replace it with the coefficient array */
    free(*offload_array);
    *offload_array = dV_drho;
    offload_data->offload_array_length = NSIZE_COMP1D * offload_data->n_rho;
 
    /* Print some sanity check on data */
    print_out(VERBOSE_IO, "\nRadial electric field (E_1DS)");
    print_out(VERBOSE_IO, "(n_rho, rho_min, rho_max) = (%d, %le, %le)\n",
              offload_data->n_rho,offload_data->rho_min,
              offload_data->rho_max);
    return err;
}
 
void E_1DS_free_offload(E_1DS_offload_data* offload_data,
                       real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void E_1DS_init(E_1DS_data* Edata, E_1DS_offload_data* offload_data,
               real* offload_array) {
    interp1Dcomp_init_spline(&Edata->dV, offload_array,
                             offload_data->n_rho,
                             NATURALBC,
                             offload_data->rho_min,
                             offload_data->rho_max);
}
 
a5err E_1DS_eval_E(real E[3], real r, real phi, real z, E_1DS_data* Edata,
                   B_field_data* Bdata) {
    a5err err = 0;
    int interperr = 0; /* If error happened during interpolation */
    real rho_drho[4];
    err = B_field_eval_rho_drho(rho_drho, r, phi, z, Bdata);
    if(err) {
        return error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_E_1DS );
    }
    /* Convert partial derivative to gradient */
    rho_drho[2] = rho_drho[2]/r;
    /* We set the field to zero if outside the profile. */
    if (rho_drho[0] < Edata->dV.x_min || rho_drho[0] > Edata->dV.x_max ) {
        E[0] = 0;
        E[1] = 0;
        E[2] = 0;
        return err;
    }
    real dV;
    interperr += interp1Dcomp_eval_f(&dV, &Edata->dV, rho_drho[0]);
 
    E[0] = -dV * rho_drho[1];
    E[1] = -dV * rho_drho[2];
    E[2] = -dV * rho_drho[3];
 
    if(interperr) {
        err = error_raise( ERR_INPUT_EVALUATION, __LINE__, EF_E_1DS );
    }
 
    return err;
}