grid_quantity.hpp

Go to the documentation of this file.

//******************************************************************************
//******************************************************************************
 
#ifndef grid_quantity_hpp
#define grid_quantity_hpp
 
#include <vector>
#include <cmath>
#include <cstddef>
 
#include "parity.hpp"
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
class radial_quantity {
public:
    const std::vector<double> buffer;
    const double ds;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    radial_quantity(const std::vector<double> &buffer);
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    virtual double get(const double s) const=0;
};
 
class half_grid;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
class full_grid : public radial_quantity {
public:
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    full_grid(const std::vector<double> &buffer);
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get(const double s) const final;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    half_grid get_prime() const;
};
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
class half_grid : public radial_quantity {
public:
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    half_grid(const std::vector<double> &buffer);
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get(const double s) const final;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    full_grid get_prime(const double m) const;
};
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <class GRID_CLASS> class vmec_grid {
public:
    const GRID_CLASS grid;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    vmec_grid(const std::vector<double> &buffer) : grid(buffer) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    vmec_grid(const GRID_CLASS &grid) : grid(grid.buffer) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get(const double s) const {
        return grid.get(s);
    }
};
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <class GRID_CLASS> class siesta_grid {
public:
    const GRID_CLASS grid;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    siesta_grid(const std::vector<double> &buffer) : grid(buffer) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    siesta_grid(const GRID_CLASS &grid) : grid(grid.buffer) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get(const double s) const {
        return grid.get(sqrt(s));
    }
};
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template<class GRID_CLASS> using vmec_quantity = std::vector<vmec_grid<GRID_CLASS> >;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <class GRID_CLASS, class PARITY> class vmec_fourier {
public:
    const std::vector<double> m;
    const std::vector<double> n;
    const vmec_quantity<GRID_CLASS> quantity;
    const PARITY func;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    vmec_fourier(const vmec_quantity<GRID_CLASS> &buffer,
                 const std::vector<double> &m,
                 const std::vector<double> &n) :
    quantity(buffer), m(m), n(n) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get(const double s,
               const double u,
               const double v) const {
        double temp = 0.0;
 
        for (size_t i = 0, e = m.size(); i < e; i++) {
            temp += quantity[i].get(s)*func.f(m[i]*u - n[i]*v);
        }
 
        return temp;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get_du(const double s,
                  const double u,
                  const double v) const {
        double temp = 0.0;
 
        for (size_t i = 0, e = m.size(); i < e; i++) {
            temp += m[i]*quantity[i].get(s)*func.df(m[i]*u - n[i]*v);
        }
 
        return temp;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get_dv(const double s,
                  const double u,
                  const double v) const {
        double temp = 0.0;
 
        for (size_t i = 0, e = m.size(); i < e; i++) {
            temp -= n[i]*quantity[i].get(s)*func.df(m[i]*u - n[i]*v);
        }
 
        return temp;
    }
};
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template<class GRID_CLASS> using siesta_quantity_1d = std::vector<siesta_grid<GRID_CLASS> >;
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template<class GRID_CLASS> using siesta_quantity = std::vector<siesta_quantity_1d<GRID_CLASS> >;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
template <class GRID_CLASS, class PARITY> class siesta_fourier {
public:
    const size_t mpol;
    const size_t ntor;
    const std::vector<int> tor_modes;
    const size_t nfp;
    const siesta_quantity<GRID_CLASS> quantity;
    const PARITY func;
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    siesta_fourier(const siesta_quantity<GRID_CLASS> &buffer,
                   const size_t mpol,
                   const size_t ntor,
                   const std::vector<int> &tor_modes,
                   const size_t nfp) :
    quantity(buffer), mpol(mpol), ntor(ntor), tor_modes(tor_modes), nfp(nfp) {}
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get(const double s,
               const double u,
               const double v) const {
        double temp = 0.0;
 
        for (size_t m = 0, em = mpol; m <= em; m++) {
            for (long n = -ntor, en = ntor; n <= en; n++) {
                const size_t ni = n + ntor;
                temp += quantity[m][ni].get(s)*func.f(m*u + tor_modes[ni]*(nfp*v));
            }
        }
 
        return temp;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get_du(const double s,
                  const double u,
                  const double v) const {
        double temp = 0.0;
 
        for (size_t m = 0, em = mpol; m <= em; m++) {
            for (long n = -ntor, en = ntor; n <= en; n++) {
                const size_t ni = n + ntor;
                temp += m*quantity[m][ni].get(s)*func.df(m*u + tor_modes[ni]*(nfp*v));
            }
        }
 
        return temp;
    }
 
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
    double get_dv(const double s,
                  const double u,
                  const double v) const {
        double temp = 0.0;
 
        for (size_t m = 0, em = mpol; m <= em; m++) {
            for (long n = -ntor, en = ntor; n <= en; n++) {
                const size_t ni = n + ntor;
                temp += tor_modes[ni]*(nfp*quantity[m][ni].get(s))*func.df(m*u + tor_modes[ni]*(nfp*v));
            }
        }
 
        return temp;
    }
};
 
#endif