BORAY: An Axisymmetric Ray Tracing Code Supports Both Closed and Open Field Lines Plasmas
(r,phi,z), axisymmetry df0/dphi=0, wave vector (kr,nphi=kphi*r,kz), cold plasma for ray tracing, hot plasma for power absorb.
Input: Br(r,z), Bz(r,z), Bphi(r,z), ns0(r,z), Ts0(r,z); and initial (r,phi,z,nphi,kz).
Support waves: Electron cyclotron wave (ECW), lower hybrid wave (LHW), ion cyclotron wave (ICW) and Alfven wave (AW).
Support configurations: Tokamak, spherical tokamak (ST), field-reversed-configuration (FRC), mirror and dipole, etc.
Cite: Hua-sheng XIE, Banerjee Debabrata, Yu-kun BAI, Han-yue ZHAO and Jing-chun LI, BORAY: An Axisymmetric Ray Tracing Code Supports Both Closed and Open Field Lines Plasmas, 2021, https://arxiv.org/abs/2105.12014.
At present, only Matlab version of BORAY is available, which has been tested in MacOS and Windows 10 under Matlab v2015b.
Files structures:
- boray_main.m % main program, run this code to start the calculation
./modules % sub files of the program
- boray_calpower.m % calculate the power absorption along the ray use hotdreach.m
- boray_initialize.m % initialize the equilibrium and other parameters
- boray_plot.m % plot the ray trajectory
- boray_plotpower.m % plot the power absorption
- boray_rk4.m % solve the cold plasma ray tracing use RK4
- calpars_solovev.m % function to calculate parameters for analytical Solovev equilibrium
- calpars.m % function to calculate parameters of numerical equilibrium
- cinterp2d.m % function to obtain bi-linear 2d interpolation coefficients
- colddr.m % to check whether the cold plasma ray trajectory is calculated accurately by solving the cold plasma dispersion relation with the k along the ray, i.e., to check whether the output omega is close to the input omega_0
- dydt.m % right hand side of the cold plasma ray tracing equation
- fcinterp.m % obtain interpolation coefficients of magnetic fields, density and temperature for numerical equilibrium
- hotdreach.m % solve the hot kinetic plasma dispersion relation to obtain omgea_i for power absorption
- initialsolovev.m % initialize the analytical Solovev equilibrium
./eqdata % input profiles data for numerical equilibrium
./output % store the output and plot figures
How to run the code:
- Choose the icase in boray_main.m.
- icase=1, with numerical EFIT equilibrium and compare with GENRAY results (store in ./eqdata directory).
- icase=2, numerical equilibrium generated by any approaches. See examples in ./eqdata directory.
- icase=3, analytical Solovev equilibrium, which supports tokamak, ST, FRC and mirror configurations.
For most cases, one can use icase=2.
-
Set initial frequency f, initial position and wave vectors yray0, which supports multi-rays.
-
Set time step dt and total # nt. Usually, larger f requires smaller dt.
-
The line ‘kr=fsolve(fDkr,kr_guess,options)’ will solve the initial kr. Use different kr_guess for different modes, e.g., X mode and O-mode.
-
Choose jcalpower for whether calculate the power absorption.
- jcalpower=0, omit the power absorption.
- jcalpower=1, calculate the power absorption.
- For power absorption part.
- jeach=0, only calculate the damping from all species.
- jeach=1, calculate the damping from each species, i.e., set other species be cold species.
-
N=4; % number of sum_N for kinetic DR, i.e., how many resonant harmonic to be kept. Increase N to make sure convergent. J=12; % =4,8,12. number of sum_J for Z(zeta) function, larger J can be accurate but also slower.
-
After seting the above parameters in boray_main.m, run it. After seconds, you can see some plots of the results. And all results will be stored in output directory. You can modify the files boray_plot.m and boray_plotpower.m for different plots, or writing your own plot files.
Of how to generate the initial numerical equilibrium profiles, see examples in ./eqdata directory. Usually, you need provide Br, Bz, Bphi, B, ns0, Ts0 in uniform (r,z) grids and also calculate the d/dr and d/dz of Br, Bz, Bphi, B, ns0.
If you meet any problems, please contact me.
Hua-sheng XIE, huashengxie@gmail.com, ENN
2021-06-02 22:52