Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
function [nel_out,s_out,nes_out,varargout] = line_average(shot,time,Rdiag,Zdiag,ne_in,rhopol_in,Rmesh,Zmesh,psi_RZmesh,psi_axis,psi_edge);
%
% compute s and ne(s) in s_out and nes_out; where s is distance from first point into the plasma along Rdiag,Zdiag chord
% compute int[ne(s)ds]/int[ds] in nel_out
%
% can get Rmesh, Zmesh and psi_RZmesh from eqdsk at shot,time from gdat if not given, but for many times, should use equil output from gdat
%
nel_out = [];
s_out = [];
nes_out = [];
time_eff = time;
if ~exist('Rmesh') || isempty(Rmesh)
% do only 1st time
time_eff = time(1);
ab=gdat(shot,'eqdsk','time',time_eff);
eqdsk=ab.eqdsk;
if isempty(Rdiag)
% assume H-1 line
par=sf2par('DCN',shot,'H1LFS_R','DCNgeo');
Rlfs=par.value/1e3;
par=sf2par('DCN',shot,'H1LFS_z','DCNgeo');
Zlfs=par.value/1e3;
par=sf2par('DCN',shot,'H1HFS_R','DCNgeo');
Rhfs=par.value/1e3;
par=sf2par('DCN',shot,'H1HFS_z','DCNgeo');
Zhfs=par.value/1e3;
Rdiag=[Rlfs Rhfs];
Zdiag=[Zlfs Zhfs];
end
if length(Rdiag)==2
nbpoints=50;
Rout=linspace(Rdiag(1),Rdiag(2),nbpoints);
Zout=linspace(Zdiag(1),Zdiag(2),nbpoints);
elseif length(Rdiag)>10
Rout=Rdiag;
Zout=Zdiag;
else
disp(['unexpected nb of diag points: length(Rdiag) = ',length(Rdiag)]);
return
end
Rmesh = eqdsk.rmesh;
Zmesh = eqdsk.zmesh;
psi_RZmesh(:,:,1) = eqdsk.psirz;
psi_axis(1) = eqdsk.psiaxis;
psi_edge(1) = eqdsk.psiedge;
end
if prod(size(Rmesh)) == max(size(Rmesh))
Rmesh_eff = repmat(reshape(Rmesh,length(Rmesh),1),1,length(time_eff));
Zmesh_eff = repmat(reshape(Zmesh,length(Zmesh),1),1,length(time_eff));
end
for it=1:length(time_eff)
psi_at_routzout(:,it) = interpos2Dcartesian(Rmesh_eff(:,it),Zmesh(:,it),psi_RZmesh(:,:,it),Rout,Zout);
ij_in=find((psi_at_routzout(:,it)-psi_axis(it)).*(psi_at_routzout(:,it)-psi_edge(it))<=0);
if ~isempty(ij_in)
if ij_in(1)==1
Rstart=Rout(ij_in(1));
Zstart=Zout(ij_in(1));
else
Rstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Rout(ij_in(1)-1:ij_in(1))) + Rout(ij_in(1)-1);
Zstart = (psi_edge(it)-psi_at_routzout(ij_in(1)-1),it)./diff(psi_at_routzout(ij_in(1)-1:ij_in(1),it)).*diff(Zout(ij_in(1)-1:ij_in(1))) + Zout(ij_in(1)-1);
end
if ij_in(end)==length(Rout)
Rend=Rout(ij_in(end));
Zend=Zout(ij_in(end));
else
Rend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Rout(ij_in(end):ij_in(end)+1)) + Rout(ij_in(end));
Zend = (psi_edge(it)-psi_at_routzout(ij_in(end),it))./diff(psi_at_routzout(ij_in(end):ij_in(end)+1,it)).*diff(Zout(ij_in(end):ij_in(end)+1)) + Zout(ij_in(end));
end
s_eff = sqrt((Rout(ij_in)-Rstart).^2 + (Zout(ij_in)-Zstart).^2);
rhopol_eff=sqrt((psi_at_routzout(ij_in,it)-psi_axis(it))./(psi_edge(it)-psi_axis(it)));
if ij_in(1)==1
rhopol_eff = reshape(rhopol_eff,length(rhopol_eff),1);
s_eff = reshape(s_eff,length(s_eff),1);
else
rhopol_eff=[1 ;reshape(rhopol_eff,length(rhopol_eff),1); 1];
s_eff=[0 ;reshape(s_eff,length(s_eff),1); sqrt((Rend-Rstart).^2 + (Zend-Zstart).^2)];
end
[rho_in_eff,isort]=sort(rhopol_in(:,it));
if rho_in_eff(1)<1e-3
% assume 1st point is center
nes = interpos(rho_in_eff,ne_in(isort,it),rhopol_eff,-0.1,[1 0],[0 0]);
else
rho_in_eff = [0 ; rho_in_eff];
ne_in_eff = [ne_in(isort(1),it); ne_in(isort,it)];
sigma_in = ones(size(ne_in_eff));
sigma_in(1)=100;
nes = interpos(rho_in_eff,ne_in_eff,rhopol_eff,-1,[1 0],[0 0],sigma_in);
end
nes_out(:,it) = nes;
s_out(:,it) = s_eff;
[dum1,~,~,nel_out_prof] = interpos(s_eff,nes,-0.1);
nel_out(it) = nel_out_prof(end)./(s_eff(end)-s_eff(1));
end
end