wall_2d.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "../ascot5.h"
#include "../print.h"
#include "wall_2d.h"
 
int wall_2d_init_offload(wall_2d_offload_data* offload_data,
                         real** offload_array) {
    // Do no initialization
 
    int n = offload_data->n;
    real rmin = (*offload_array)[0], rmax = (*offload_array)[0];
    real zmin = (*offload_array)[n], zmax = (*offload_array)[n];
    for(int i=0; i<n; i++) {
        rmin = fmin(rmin, (*offload_array)[i]);
        rmax = fmax(rmax, (*offload_array)[i]);
        zmin = fmin(zmin, (*offload_array)[n+i]);
        zmax = fmax(zmax, (*offload_array)[n+i]);
    }
 
    print_out(VERBOSE_IO, "\n2D wall model (wall_2D)\n");
    print_out(VERBOSE_IO, "Number of wall elements = %d,"
              " R extend = [%2.2f, %2.2f], z extend = [%2.2f, %2.2f]\n",
              n, rmin, rmax, zmin, zmax);
 
    return 0;
}
 
void wall_2d_free_offload(wall_2d_offload_data* offload_data,
                          real** offload_array) {
    free(*offload_array);
    *offload_array = NULL;
}
 
void wall_2d_init(wall_2d_data* w, wall_2d_offload_data* offload_data,
                  real* offload_array) {
    w->n = offload_data->n;
    w->wall_r = &offload_array[0];
    w->wall_z = &offload_array[offload_data->n];
}
 
int wall_2d_inside(real r, real z, wall_2d_data* w) {
    int hits = 0;
    for(int i = 0; i < w->n; i++) {
        real wr1, wr2, wz1, wz2;
        if(i == w->n - 1) {
            wz1 = w->wall_z[i] - z;
            wz2 = w->wall_z[0]   - z;
            wr1 = w->wall_r[i] - r;
            wr2 = w->wall_r[0]   - r;
        } else {
            wz1 = w->wall_z[i]   - z;
            wz2 = w->wall_z[i+1] - z;
            wr1 = w->wall_r[i]   - r;
            wr2 = w->wall_r[i+1] - r;
        }
        if(wz1 * wz2 < 0) {
            real ri = wr1 + (wz1*(wr2-wr1)) / (wz1-wz2);
            if(ri > 0) {
                hits++;
            }
        }
    }
    return hits % 2;
}
 
int wall_2d_hit_wall(real r1, real phi1, real z1, real r2, real phi2, real z2,
                     wall_2d_data* w, real* w_coll) {
    int tile = 0;
    if(!wall_2d_inside(r2, z2, w)) {
        tile = wall_2d_find_intersection(r1, z1, r2, z2, w, w_coll);
    }
    return tile;
}
 
int wall_2d_find_intersection(real r1, real z1, real r2, real z2,
                              wall_2d_data* w, real* w_coll) {
    int tile = 0;
    real t0 = 2.0; // Helper variable to pick the closest intersection
    for(int i=0; i<w->n; i++) {
        real r3, z3, r4, z4;
        if(i == w->n-1) {
            r3 = w->wall_r[i];
            z3 = w->wall_z[i];
            r4 = w->wall_r[0];
            z4 = w->wall_z[0];
        } else {
            r3 = w->wall_r[i];
            z3 = w->wall_z[i];
            r4 = w->wall_r[i+1];
            z4 = w->wall_z[i+1];
        }
 
        real div = (r1 - r2) * (z3 - z4) - (z1 - z2) * (r3 - r4);
        real t   = ( (r1 - r3) * (z3 - z4) - (z1 - z3) * (r3 - r4) ) / div;
        real u   = ( (r1 - r3) * (z1 - z2) - (z1 - z3) * (r1 - r2) ) / div;
        if(0 <= t && t <= 1.0 && 0 <= u && u <= 1.0 && t < t0) {
            t0   = t;
            tile = i + 1;
        }
    }
    *w_coll = t0;
    return tile;
}