Updates

function [alpha_deg,beta_deg] = calc_Bangle_g(g,R1,Z1,R2,Z2)

error('verify that midpoint is desired method')

%  alpha_deg = calc_Bangle_g(g,R1,Z1,R2,Z2)

function P = calc_Bangle_g(g,RZ1,RZ2,npts)

%

% RZ1, RZ2 ... start and end points used to calculate surface normal.

% Surface normal is (RZ2-RZ1,0)x(0,0,1) for (R,Z,phi), with phi coming out

% of the page. This means if wall elements are anticlockwise then the

% surface normal will point "out", away from the axis.

%

% npts is number of points to evaluate along this segment. 

% If npts == 1 the midpoint is used.

%

% P.alpha is the angle of incidence relative to the surface, i.e., the 

% incident flux is qdep = qprl*sin(P.alpha)

%

% P.beta is the angle in the poloidal plane between B and the surface

% normal

%

if length(RZ1) ~= 2 || length(RZ2) ~= 2

    error('Only single point pairs are allowed')

end

if nargin < 3

    npts = 3;

end

if npts < 1

    error('npts must be at least 1')

end

R1 = RZ1(1);

Z1 = RZ1(2);

R2 = RZ2(1);

Z2 = RZ2(2);

DEBUG = 0;

if DEBUG

    dPlot = 0.01;

    figure; hold on; box on; set(gcf,'color','w');set(gca,'fontsize',14,'fontweight','bold'); grid on;

    axis equal;

    plot([R1,R2],[Z1,Z2],'ko-')

end

v2 = [0,0,1];

for i = 1:length(R1)

    v1 = [R2(i)-R1(i),Z2(i)-Z1(i),0];

if npts == 1

    Reval = 0.5*(R2+R1);

    Zeval = 0.5*(Z2+Z1);

else

    Reval = linspace(R1,R2,npts);

    Zeval = linspace(Z1,Z2,npts);

end

if DEBUG 

    plot(Reval,Zeval,'rx')

end

% Surface normal is the same across the segment

v1 = [Reval(end)-R1,Zeval(end)-Z1,0];

vn = cross(v1,v2);

vn = vn./norm(vn);

    b = bfield_geq_bicub(g,0.5*(R2(i)+R1(i)),0.5*(Z2(i)+Z1(i)),0);

% Then local B vector is used

alpha_deg = zeros(1,npts);

beta_deg = zeros(1,npts);

for j = 1:length(Reval)

    if DEBUG

        t = atan2(vn(2),vn(1));

        plot([Reval(j),Reval(j)+dPlot*cos(t)],[Zeval(j),Zeval(j)+dPlot*sin(t)],'b')        

    end

    b = bfield_geq_bicub(g,Reval(j),Zeval(j),0);

    b = [b.br,b.bz,b.bphi];

    b = b./norm(b);

    alpha_deg(i) = 90 - acos(dot(vn,b))*180/pi;

    if DEBUG

        t = atan2(b(2),b(1));

        plot([Reval(j),Reval(j)+dPlot*cos(t)],[Zeval(j),Zeval(j)+dPlot*sin(t)],'m')        

    end

    alpha_deg(j) = asin(dot(vn,b))*180/pi;

    brz = b;

    brz(3) = 0;

    brz = brz./norm(brz);

    beta_deg(i) = acos(dot(vn,brz))*180/pi;

    beta_deg(j) = acos(dot(vn,brz))*180/pi;

end

end

 No newline at end of file

P.alpha_deg = alpha_deg;

P.beta_deg = beta_deg;

function sep = plot_sep_g(g,plotit,newfig,linewid)

function sep = plot_sep_g(g,plotit,newfig,linewid,clip_at_glim,flush_file)

% function sep = plot_sep_g(g,plotit,newfig)

if nargin < 2

    plotit = 1;

if nargin < 4

    linewid = 1;

end

if nargin < 5

    clip_at_glim = 1;

end

if nargin < 6

    flush_file = 0;

end

if newfig == 1

    figure; hold on;   

end

DEBUG = 0;

box_r = [g.r(1),g.r(g.mw),g.r(g.mw),g.r(1),g.r(1)];

[pathstr,fn,ext] = fileparts(g.filename);

fname = fullfile(pathstr,strcat(fn,ext,'_sep'));

if flush_file && isfile(fname)

    delete(fname);

end

if DEBUG 

    figure; hold on;

    Z1 = Z1.';

end

% This is a switch to return B = Bt = 1 when the evaluation point(s)

% are off the grid. 

toroidal_off_grid = 0;

if isfield(g,'toroidal_off_grid')

    toroidal_off_grid = g.toroidal_off_grid;

clearvars;

gfile_name = 'C:\Users\jjl\Dropbox (ORNL)\DIII-D\Qprl experiment\power13mw\g174308.03500_159';

g = readg_g3d(gfile_name);

Rtest = linspace(1.8,3.4,100);

Ztest = linspace(-1,1,100);

bb = bfield_geq_bicub_inv(g,Rtest,Ztest);

% fprintf('R, Z, Br, Bphi, Bz = [%8.6f, %8.6f, %8.6f, %8.6f, %8.6f]\n',Rtest,Ztest,bb.br,bb.bphi,bb.bz)