Skip to content
Snippets Groups Projects
Commit 7ec5b296 authored by Olivier Sauter's avatar Olivier Sauter
Browse files

add files needed for sxr and ECE in TCV/private mainly 

git-svn-id: https://spcsvn.epfl.ch/repos/TCV/gdat/trunk@1823 d63d8f72-b253-0410-a779-e742ad2e26cf
parent bbf70925
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
TCV/VsxrTCVradius.m 0 → 100644
+ 24
− 0
View file @ 7ec5b296
function radius=VsxrTCVradius(ymag,xchord,ychord)
% get intersection of magnetics axis with chords
% input:
% Vsawtooth.shot : shot number
% Vsawtooth.expdata.zmag : positions of the magnetic a
% xchord: two x-coordinates.
% ychord: two y-coordinates.
% For each line (2xnl), they specify start + end pointsxis
% Output :
% Vsxr.radius : intersectionof magnetics axis with chords
%parameter: start and end of each of the camera number 2
% calculation intersections
xchord=xchord/100;
ychord=ychord/100;
for i=1:size(xchord,2)
a=(ychord(1,i)-ychord(2,i))/(xchord(1,i)-xchord(2,i));
b=ychord(2,i)-a*xchord(2,i);
radius(:,i)= (ymag-b)/(a+eps);
end
TCV/ece_te.m 0 → 100755
+ 67
− 0
View file @ 7ec5b296
function [TE_ECE,TE_ECE_ERR,RHO,R,T,TE_THOM,TE_THOM_ERR,Fcentral,CAL,NE_THOM]=ece_te(shot,Tc,T,SS)
% [TE_ECE,TE_ECE_ERR,RHO,R,T,TE_THOM,TE_THOM_ERR,Fcentral,CAL,NE_THOM]=ece_te(shot,Tc,T,SS)
% [TE_ECE,TE_ECE_ERR,RHO,R,T,TE_THOM,TE_THOM_ERR,Fcentral,CAL,NE_THOM]=ece_te(19314,[0.1 0.29],[0 1.5],10);
%
% Program that compute TE_ECE and TE_THOM profile over R or RHO
% for all the times T with the error matrices TE_ECE_ERR and TE_THOM_ERR.
% The calibration matrice CAL comes from the ece_calib.m program.
% Fcentral is the central frequency of the ECE working channels.
%
% shot = shot number
% Tc = [a b] temporal boundary on which we want the calibration on
% thomson temperature. If you have no idea, put Tc=[0.1 0.29]
% T = either
% -temporal vector
% -if T=10 => takes all the ECE times
% -if T=[c d];=> c and d are the temporal boundary for the output
% SS >= 1 => sampling rate SS if T=10 or T=[c d]
%
% B est calcule avec la routine BandBres_allt de O.Sauter
%
% Blanchard 25.11.2000
%-----------------------------------------------------------------------
% Recherche de la configuration de l'ECE
%-----------------------------------------------------------------------
[RHOece,Rece,Zece,Tece,Fcentralrho]=ece_rho(shot,T);
[CAL,Fcentral]=ece_calib(shot,Tc);
%-----------------------------------------------------------------------
% Recherches des signaux voulus
%-----------------------------------------------------------------------
[ECE,TECE,Fcentral]=ece_raw_signals(shot,T,SS);
%[RHOece,Rece,Zece,Tece, Fcentralrho]=ece_rho(shot,T);
%Zece=Zece*ones(size(Rece));
[TEthom,NEthom,TEerr,NEerr,Tthom,RHOthom]=thom_rho(shot,6);
if exist('TEthom')==0|length(TEthom) ==0
disp('Le profil Thomson proffit n''existe pas. On prend le profil direct') %(Tthom,Nbre de pts selon rho)
[TEthom,NEthom,TEerr,NEerr,Tthom,RHOthom]=thom_rho(shot,3); %(Tthom,Nbre de pts selon rho)
end
[t2,i3,i4]=common_ece(Fcentralrho,Fcentral);
Rece=Rece(i3,:);Zece=Zece(i3,:);RHOece=RHOece(i3,:);
TE_ECE=ECE.*repmat(CAL(2,:),length(TECE),1);
T=repmat(TECE,1,length(Fcentral));
if length(Tece)==1
RHO=repmat(RHOece,1,length(TECE));RHO=RHO';
R=repmat(Rece,1,length(TECE));R=R';
X=size(RHOthom);
else
%RHO=interp2(repmat(Fcentral,1,length(Tece))',repmat(Tece,length(Fcentral),1)', ...
% RHOece',repmat(Fcentral',length(TECE),1)',T');
RHO=interp2(repmat(Fcentral,1,length(Tece))',repmat(Tece',length(Fcentral),1)', ...
RHOece',repmat(Fcentral',length(TECE),1)',T');
RHO=RHO';
R=interp2(repmat(Fcentral,1,length(Tece))',repmat(Tece',length(Fcentral),1)', ...
Rece',repmat(Fcentral',length(TECE),1),T);
end
X=size(RHOthom);
TE_THOM=griddata(repmat(Tthom,1,X(2)),RHOthom,TEthom,T,RHO);
NE_THOM=griddata(repmat(Tthom,1,X(2)),RHOthom,NEthom,T,RHO);
TE_THOM_ERR=griddata(repmat(Tthom,1,X(2)),RHOthom,TEerr,T,RHO);
TE_ECE_ERR=CAL(4,:);T=T(:,1);
TCV/get_xtomo_data.m 0 → 100755
+ 497
− 0
View file @ 7ec5b296
function [sig,t]=get_xtomo_data(shot,t1,t2,dt,fans,angfact,tag);
% -
%[sig,t]=get_xtomo_data_m5(shot,t1,t2,dt,fans,angfact,tag);
%
% INPUT:
% shot: TCV shot
% t1: start time
% t2: stop time
% dt: timestep
% fans: camera switch, e.g. [0 0 0 0 0 0 1 0 1 0];
% angfact: relative etendue, size: [20 x 10]
% tag 'full' or 'brief', indicates if it is
% a slow or a fast Pentland acquisition
%
% OUTPUTS:
% sig: xtomo signals, size: [sum(fans) x length(t)]
% t: times
%
% ATTENTION: length(time) may be shorter than foreseen !!
%
% This routine works on Matlab5.
% Original routine for Matlab4 by Anton Mathias.
%
% Last update: 25-08-1999
%
%-------------MAC:[FURNO.MATLAB5.XTOMO]----------------------------------
%---- get data and offsets ---------------------------------------------------
% set a flag if flattenign of noise
% for channels with low signal
iflat=0;
minsiglevel=0.01;
satlevel=9.9;
st1=sprintf('%6.4f',t1);
st2=sprintf('%6.4f',t2);
sdt=num2str(dt);
% sdt=sprintf('%6.4f',dt); % modified from old version
if shot >= 13836 & shot <= 13848
tstart=-0.02;tstop=-0.01; % this one is to be used only for
else % shot=13836 to shot=13848
tstart=-0.04;tstop=-0.01;
end
ststart=num2str(tstart);
ststop=num2str(tstop);
if nargin <=6
shot=mdsopen('eltca1::tcv_shot',shot);
t=mdsdata(['dim_of(\base::xtomo:array_001[',st1,':',st2,':',sdt,',*])']);
S1=[];
S2=[];
S3=[];
S4=[];
S5=[];
S6=[];
S7=[];
S8=[];
S9=[];
S10=[];
if isempty(t)
disp('get_xtomo_data: sorry, nothing to be found for this shot ...')
return
else
disp('*-------------------------------------------*')
disp('| get_xtomo_data: getting data from MDS |')
disp('*-------------------------------------------*')
end
if shot>6768
if fans(1)
S=mdsdata(['\base::xtomo:array_001[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_001[',ststart,':',ststop,',0:19]']);
S1=S-repmat(mean(offset),length(t),1)';
end
if fans(2)
S=mdsdata(['\base::xtomo:array_002[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_002[',ststart,':',ststop,',0:19]']);
S2=S-repmat(mean(offset),length(t),1)';
end
if fans(3)
S=mdsdata(['\base::xtomo:array_003[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_003[',ststart,':',ststop,',0:19]']);
S3=S-repmat(mean(offset),length(t),1)';
end
if fans(4)
S=mdsdata(['\base::xtomo:array_004[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_004[',ststart,':',ststop,',0:19]']);
S4=S-repmat(mean(offset),length(t),1)';
end
if fans(5)
S=mdsdata(['\base::xtomo:array_005[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_005[',ststart,':',ststop,',0:19]']);
S5=S-repmat(mean(offset),length(t),1)';
end
if fans(6)
S=mdsdata(['\base::xtomo:array_006[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_006[',ststart,':',ststop,',0:19]']);
S6=S-repmat(mean(offset),length(t),1)';
end
if fans(7)
S=mdsdata(['\base::xtomo:array_007[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_007[',ststart,':',ststop,',0:19]']);
S7=S-repmat(mean(offset),length(t),1)';
end
if fans(8)
S=mdsdata(['\base::xtomo:array_008[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_008[',ststart,':',ststop,',0:19]']);
S8=S-repmat(mean(offset),length(t),1)';
end
if fans(9)
S=mdsdata(['\base::xtomo:array_009[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_009[',ststart,':',ststop,',0:19]']);
S9=S-repmat(mean(offset),length(t),1)';
end
if fans(10)
S=mdsdata(['\base::xtomo:array_010[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_010[',ststart,':',ststop,',0:19]']);
S10=S-repmat(mean(offset),length(t),1)';
end
sig=-[S1;S2;S3;S4;S5;S6;S7;S8;S9;S10];
[satrow,satcol]=find(abs(sig)>satlevel);
if ~isempty(satcol)
i_tlimit=min(satcol(:));
if i_tlimit>1
sig=sig(:,1:i_tlimit-1);
t=t(1:i_tlimit-1);
disp(['get_xtomo_data WARNING: some channels saturated',......
', t2 changed to ',sprintf('%6.4f',max(t))]);
else
sig=[];
t=[];
disp('get_xtomo_data WARNING: saturations, no data returned');
return
end
end
angfact=angfact(:,find(fans));
angfact=angfact(:);
gains=get_xtomo_gains(shot);
% earlier than shot xxx:
perm=[1:180,[180+[2,1,4,3,6,5,8,7,10,9,12,11,14,13,16,15,18,17,20,19]]];
gains=gains(perm);
iact=find(fans);
nact=length(iact);
idetec=[];
for k=1:nact
idetec=[idetec,(iact(k)-1)*20+1:iact(k)*20];
end
gains=gains(idetec);
gains=gains(:);
fac=ones(size(gains))./gains;
sig=-sig.*(repmat(fac,1,length(t)));
if iflat
for ii=1:length(t)
imini=find( sig(:,ii) < minsiglevel*max(sig(:,ii)));
sig(imini,ii)=zeros(size(imini));
end
end
t=t';
sig=sig.*(repmat(angfact,1,length(t)));
else
if fans(9)
S=mdsdata(['\base::xtomo:array_001[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_001[',ststart,':',ststop,',0:19]']);
V=S-repmat(mean(offset),length(t),1)';
end
if fans(7)
S=mdsdata(['\base::xtomo:array_002[',st1,':',st2,':',sdt,',0:19]'])';
offset=mdsdata(['\base::xtomo:array_002[',ststart,':',ststop,',0:19]']);
H=S-repmat(mean(offset),length(t),1)';
end
%---- read gains of the current shot ---------------------------------------
gains=get_xtomo_gains(shot);
g1=gains(1:20)';
g2=gains(21:40)';
%------ solid angle factors -----------------------
aom1=angfact(:,9);
aom2=angfact(:,7);
%----- calculate correcting factors ---------------------------------------
fac1=(aom1./g1);
fac2=(aom2./g2);
%--- result: offset-, calibration- and angular-factor corrected data -------
V=-V.*(repmat(fac1,1,length(t))); % vertical camera:
H=-H.*(repmat(fac2,1,length(t))); % horizontal camera:
sig=[H;V];
t=t';
end
else
if (~strcmp(tag,'brief') & ~strcmp(tag,'full'))
disp('Only full or brief accepted as tag')
return
end
trace_tree=['\atlas::t_rex3_' tag ':'];
trace=[trace_tree 'signal_inp:p'];
%shot=mdsopen('eltca1::tcv_shot',shot);
shot=mdsopen(shot);
tref=mdsdata(['_tref=dim_of(' trace,'001)']); %trace :p001
dtref=mdsdata('_dtref = _tref[1]-_tref[0]');
S1=[];
S2=[];
S3=[];
S4=[];
S5=[];
S6=[];
S7=[];
S8=[];
S9=[];
S10=[];
t_1=mdsdata('_t = dim_of(\atlas::t_rex3_full:signal_inp:p001)');
dt1=mdsdata('_dt = _t[1]-_t[0]');
npts=floor(mdsdata(['_npts=(',ststop,'-(',ststart,'))/_dt']));
if isempty(t_1)
disp('get_xtomo_data: sorry, nothing to be found for this shot ...')
return
else
disp('*--------------------------------------------*')
disp('| getting Pentland acquisition: getting *')
if strcmp(tag,'brief')
disp('| BRIEF data from MDS *');
else
disp('| FULL data from MDS *');
end
disp('*--------------------------------------------*');
end
if ((dt - dtref) < 1e-7 | floor(dt/dtref) ==0)
step='1';
else
step =num2str(dt/dtref);
end
ind1=find( abs(tref-t1) == min(abs(tref-t1)));
npts1=floor(mdsdata(['_npts1=(',st2,'-(',st1,'))/_dtref']));
sind1=num2str(ind1);
sind2=num2str(npts1+ind1);
t=mdsdata(['data(_tref)[',sind1,':',sind2,':', step, ']']);
if fans(1)
channel=['001';'002';'003';'004';'005';'006';'007';'008';'009';'010';...
'011';'012';'013';'014';'015';'016';'017';'018';'019';'020'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':' , step,']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S1(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(2)
channel=['021';'022';'023';'024';'025';'026';'027';'028';'029';'030';...
'031';'032';'033';'034';'035';'036';'037';'038';'039';'040'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step, ']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S2(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(3)
channel=['041';'042';'043';'044';'045';'046';'047';'048';'049';'050';...
'051';'052';'053';'054';'055';'056';'057';'058';'059';'060'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step, ']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S3(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(4)
channel=['061';'062';'063';'064';'065';'066';'067';'068';'069';'070';...
'071';'072';'073';'074';'075';'076';'077';'078';'079';'080'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step, ']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S4(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(5)
channel=['081';'082';'083';'084';'085';'086';'087';'088';'089';'090';...
'091';'092';'093';'094';'095';'096';'097';'098';'099';'100'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':',step,']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S5(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(6)
channel=['101';'102';'103';'104';'105';'106';'107';'108';'109';'110';...
'111';'112';'113';'114';'115';'116';'117';'118';'119';'120'];
for i=1:20
% S(i,:)=mdsdata([trace,channel(i,:),'[',st1,':',st2,':',sdt,':',step,']'])'
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step, ']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S6(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(7)
channel=['121';'122';'123';'124';'125';'126';'127';'128';'129';'130';...
'131';'132';'133';'134';'135';'136';'137';'138';'139';'140'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':',step,']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S7(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(8)
channel=['141';'142';'143';'144';'145';'146';'147';'148';'149';'150';...
'151';'152';'153';'154';'155';'156';'157';'158';'159';'160'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step,']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S8(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(9)
channel=['161';'162';'163';'164';'165';'166';'167';'168';'169';'170';...
'171';'172';'173';'174';'175';'176';'177';'178';'179';'180'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step,']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S9(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
if fans(10)
% channel=['181';'182';'183';'184';'185';'186';'187';'188';'189';'190';...
% '191';'192';'193';'194';'195';'196';'197';'198';'199';'200'];
% this permutation takes into account some not well defined (not yet)
% hardware corrections. Are the gains to be permutated ?
channel=['182';'181';'184';'183';'186';'185';'188';'187';'190';'189';...
'192';'191';'194';'193';'196';'195';'198';'197';'200';'199'];
for i=1:20
% S(i,:)=mdsdata([trace, channel(i,:),'[',st1,':',st2,':',sdt,']'])';
S(i,:)=mdsdata(['data(',trace,channel(i,:),')[',sind1,':',sind2,':', step,']'])';
offset=mdsdata(['data(\atlas::t_rex3_full:signal_inp:p',channel(i,:),')[0:',num2str(npts),']' ]);
S10(i,:)=S(i,:)-repmat(mean(offset),length(t),1)';
end
end
sig=[S1;S2;S3;S4;S5;S6;S7;S8;S9;S10];
[satrow,satcol]=find(abs(sig)>satlevel);
if ~isempty(satcol)
i_tlimit=min(satcol(:));
if i_tlimit>1
sig=sig(:,1:i_tlimit-1);
t=t(1:i_tlimit-1);
disp(['get_xtomo_data WARNING: some channels saturated',......
', t2 changed to ',sprintf('%6.4f',max(t))]);
else
sig=[];
t=[];
disp('get_xtomo_data WARNING: saturations, no data returned');
return
end
end
angfact=angfact(:,find(fans));
angfact=angfact(:);
gains=get_xtomo_gains(shot);
% earlier than shot xxx:
perm=[1:180,[180+[2,1,4,3,6,5,8,7,10,9,12,11,14,13,16,15,18,17,20,19]]];
gains=gains(perm);
iact=find(fans);
nact=length(iact);
idetec=[];
for k=1:nact
idetec=[idetec,(iact(k)-1)*20+1:iact(k)*20];
end
gains=gains(idetec);
gains=gains(:);
fac=ones(size(gains))./gains;
sig=sig.*(repmat(fac,1,length(t)));
if iflat
for ii=1:length(t)
imini=find( sig(:,ii) < minsiglevel*max(sig(:,ii)));
sig(imini,ii)=zeros(size(imini));
end
end
t=t';
sig=sig.*(repmat(angfact,1,length(t)));
end
mdsclose
return
TCV/private/angular_fact_1.mat 0 → 100755
+ 0
− 0
View file @ 7ec5b296
File added
TCV/private/bad_channels.mat 0 → 100755
+ 0
− 0
View file @ 7ec5b296
File added
TCV/private/ece_d.m 0 → 100755
+ 51
− 0
View file @ 7ec5b296
% Fonction matlab donnant la relation entre la frequence centrale des
% cannaux du radiometre avec le no du cannal
D=[ 78.8500 12.0000
79.6000 13.0000
80.3500 11.0000
81.1000 14.0000
81.8500 10.0000
82.6000 15.0000
83.3500 9.0000
84.1000 16.0000
84.8500 8.0000
85.6000 17.0000
86.3500 7.0000
87.1000 18.0000
87.8500 6.0000
88.6000 19.0000
89.3500 5.0000
90.1000 20.0000
90.8500 4.0000
91.6000 21.0000
92.3500 3.0000
93.1000 22.0000
93.8500 2.0000
94.6000 23.0000
95.3500 1.0000
96.1000 24.0000
96.8500 12.0000
97.6000 13.0000
98.3500 11.0000
99.1000 14.0000
99.8500 10.0000
100.6000 15.0000
101.3500 9.0000
102.1000 16.0000
102.8500 8.0000
103.6000 17.0000
104.3500 7.0000
105.1000 18.0000
105.8500 6.0000
106.6000 19.0000
107.3500 5.0000
108.1000 20.0000
108.8500 4.0000
109.6000 21.0000
110.3500 3.0000
111.1000 22.0000
111.8500 2.0000
112.6000 23.0000
113.3500 1.0000
114.1000 24.0000 ];
TCV/private/tcv_vesc1.mat 0 → 100644
+ 0
− 0
View file @ 7ec5b296
File added
TCV/xtomo_geometry.m 0 → 100644
+ 402
− 0
View file @ 7ec5b296
function [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans)
% ----[anton.public]
%
% function
%
% [fans,vangle,xchord,ychord,aomega,angfact]=xtomo_geometry(i_detec,fans);
%
% inputs:
%
% i_detec: =2: Xtomo prototype cameras (shot# < 6768)
% =1: Xtomo 9-cameras (shot# > 682x)
%
% outputs:
%
% fans: camera switch, 1=on,0=off (1x10)
% vangle: angle between detect. surface normal and pos. x-axis (1x10)
% xchord: two x-coordinates (2xnl) and
% ychord: two y-coord. for each line (2xnl), they specify start + end points
% aomega: etendue in mm^2 x steradians
% angfact: angular factors, inverse of relative etendue (throughput) (20x10)
%
% uses:
% AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw);
% angular_fact_*.mat , '*'=i_detec
%
%
%---------------- M.Anton 14/3/95 -------------------------------------------
disp('*----------------------------*')
disp('| this is xtomo_geometry |')
disp('*----------------------------*')
global xap yap xdet ydet
global ae da
% ======== tokamak parameters ================================================
load tcv_vesc1
xcont=rzvin(:,1);
ycont=rzvin(:,2);
xmin=min(xcont);
ymin=min(ycont);
xmax=max(xcont);
ymax=max(ycont);
xedge=100;
yedge=60;
% ========= detector parameters =============================================
if i_detec==2
cw=1; % detector numbers cw=1:clockwise cw=0:ccw
if nargin<2
fans=[0 0 0 0 0 0 1 0 1 0]; % camera switch
end
vangle=[90 90 90 0 0 0 0 -90 -90 -90];
% angle of detector surface normal
xpos=[0 0 0 0 0 0 118.05 0 87.84 0];
% x position of the diaphragmas in [cm]
ypos=[0 0 0 0 0 0 -46 0 -80.45 0];
% y position of the diaphragmas in [cm]
ae=[0 0 0 0 0 0 -2.5 0 -0.1 0]/10;
%excentricity array/diaphragm in [cm]
da=[0 0 0 0 0 0 10.1 0 -11.7 0]/10;
% diaphragma-array distance in [cm]
da2=da;
d1=0.950; % detector width in mm
d2=4.050; % detector length in mm
b1=1.000; % aperture width in mm
b2=4.000*ones(1,10); % aperture length in mm
b3x=0; % aperture thickness in mm
b3y=0;
elseif i_detec==1
cw=1; % detector numbers cw=1:clockwise cw=0:ccw
if nargin<2
fans=[1 1 1 1 1 1 1 1 0 1]; % camera switch
end
vangle=[90 90 90 0 0 0 0 -90 -90 -90];
% angle of detector surface normal
xpos=[71.5 87.7 103.9 123.1 123.1 123.1 123.1 104.04 87.84 71.64];
% x position of the diaphragmas in [cm]
ypos=[80.35 80.35 80.35 48.1 1.95 -2.45 -48.6 -80.35 -80.35 -80.35];
% y position of the diaphragmas in [cm]
ae=[-8 0 8 5 9 -9 -5 8 0 -8]/10; %excentricity array/diaphragm [cm]
ae= ae + [-0.0915 0 0.1361 0.2123 0.0923 -0.0994 ...
0.0872 -0.1520 0 0.9410 ]/10;
ae(1)=ae(1)+0.1/10;
ae(3)=8/10+0.14/10;
ae(4)=4.9/10;
ae(5)=9/10+0.2/10;
ae(6)=ae(6)-0.2/10;
ae(7)=-4.9/10;
ae(10)=-7.1/10;
da= [12.4 9.9 12.4 9.9 13.4 13.4 9.9 -12.4 -9.9 -12.4]/10;
% diaphragma-array distance in [cm] (poloidal)
da2=[37 34.4 37 55.9 59.4 59.4 55.9 37 34.4 37]/10;
% dist to diaphragm in toroidal direction [cm];
deltada=[ -0.0311 0 -0.0458 -0.1179 -0.0615 -0.1105 ...
-0.0510 -0.0515 0 -0.3223]/10;
deltada(4)=0;
deltada(6)=0;
da=da+ deltada;
da2=da2+deltada;
d1=0.90; % detector width in mm
d2=4.0; % detector length in mm
b1=0.800; % aperture width in mm (pol.)
b2=[8 8 8 15 15 15 15 8 8 8];
% aperture length in mm (tor.)
b3x=0.020; % aperture thickness in mm (poloidal)
b3y=0; % aperture thickness in mm (toroidal)
end
%======== calculation of the chords of view ===================================
nact=sum(fans);
iact=find(fans);
ndet=20;
ncam=10;
% ---- apertures: ------------------
xap=ones(ndet,1)*xpos(iact);
xap=xap(:)';
yap=ones(ndet,1)*ypos(iact);
yap=yap(:)';
% ---- detectors: ------------------
vorz(find(vangle==90))=(-1)^(cw+1)*ones(size(find(vangle==90)));
vorz(find(vangle==0))=(-1)^cw*ones(size(find(vangle==0)));
vorz(find(vangle==-90))=(-1)^cw*ones(size(find(vangle==-90)));
dete=(-9.025:0.950:9.025)'/10*vorz(iact)+ones(ndet,1)*ae(iact);
dum_ae=dete(:)';
dum_vangle=ones(ndet,1)*vangle(iact);
dum_vangle=dum_vangle(:)';
ivert=find(dum_vangle==90 | dum_vangle==-90);
ihori=find(dum_vangle==0);
dum_da=ones(ndet,1)*da(iact);
dum_da=dum_da(:)';
dxd=zeros(1,ndet*nact);
dyd=zeros(1,ndet*nact);
dxd(ivert)=dum_ae(ivert);
dxd(ihori)=dum_da(ihori);
dyd(ivert)=dum_da(ivert);
dyd(ihori)=dum_ae(ihori);
xdet=xap+dxd;
ydet=yap+dyd;
%plot_vessel(rzvin,rzvout)
%hold on
% plot(xap,yap,'.g',xdet,ydet,'.m')
% ---- calculate the equations of lines of sight
m=(ydet-yap)./(xdet-xap);
b=ydet-m.*xdet;
nl=length(xdet);
xchord=zeros(2,nl);
ychord=zeros(2,nl);
xchord(1,:)=xdet;ychord(1,:)=ydet;
iup=find(dum_vangle==90);
isi=find(dum_vangle==0);
ido=find(dum_vangle==-90);
if ~isempty(iup)
ychord(2,iup)=ymin*ones(size(iup));
xchord(2,iup)=(ychord(2,iup)-b(iup))./m(iup);
end
if ~isempty(ido)
ychord(2,ido)=ymax*ones(size(ido));
xchord(2,ido)=(ychord(2,ido)-b(ido))./m(ido);
end
if ~isempty(isi)
xchord(2,isi)=xmin*ones(size(isi));
ychord(2,isi)=m(isi).*xchord(2,isi)+b(isi);
end
ileft=find(xchord(2,:)<xmin);
if ~isempty(ileft)
xchord(2,ileft)=xmin*ones(size(ileft));
ychord(2,ileft)=m(ileft).*xchord(2,ileft)+b(ileft);
end
irig=find(xchord(2,:)>xmax);
if ~isempty(irig)
xchord(2,irig)=xmax*ones(size(irig));
ychord(2,irig)=m(irig).*xchord(2,irig)+b(irig);
end
%======== prepare output ======================================================
vangle=vangle(iact);
%======== calculation of angular correction factors, if necessary =============
if i_detec==2 & exist('angular_fact_2.mat')==2
disp('loading angular_fact_2')
load angular_fact_2
elseif i_detec==1 & exist('angular_fact_1.mat')==2
disp('loading angular_fact_1')
load angular_fact_1
else
aomega=zeros(ndet,ncam);
angfact=ones(ndet,ncam);
for l=1:sum(fans)
% Z0X=abs(da(iact(l))*10)
% Z0Y=abs(da2(iact(l))*10)
% X0=ae(iact(l))*10 % back to mm, sorry about that...
% X0=X0*vorz(iact(l))
% B2=b2(iact(l))
% AOMEGA=etendue_n(b1,B2,b3x,b3y,Z0X,Z0Y,X0,cw);
b1x=0.8;
b1y=6;
b1z=0.02;
b2x=10000000;
b2y=b2(iact(l));
b2z=0;
z01=abs(da(iact(l))*10);
z02=abs(da2(iact(l))*10);
X0=ae(iact(l))*10*vorz(iact(l));
keyboard
AOMEGA=etendue_n2(b1x,b1y,b1z,b2x,b2y,b2z,z01,z02,X0,cw);
aomega(:,iact(l))=AOMEGA(:,1);
end
indm=min(find(aomega==max(aomega(:))));
aomegan=aomega/aomega(indm);
nonz=find(aomega);
angfact(nonz)=ones(size(nonz))./aomegan(nonz);
angfact=round(1000*angfact)/1000;
unitstring='units aomega: mm^2 * sterad';
if i_detec==1
save angular_fact_1 angfact aomega unitstring
elseif i_detec==2
save angular_fact_2 angfact aomega unitstring
end
end
return
th=atan(diff(ychord)./diff(xchord));
thp=th;
neg=find(thp<0);
thp(neg)=180+thp(neg);
thdet=ones(1,20*sum(fans));
for k=1:sum(fans)
thdet((k-1)*20+1:k*20)=vangle(k)*ones(1,20);
end
angles=thdet-thp;
mist=find(angles<0 & abs(angles)>90);
angles(mist)=angles(mist)+180;
th_inc=angles*pi/180;
% ---- correct for the edges of tcv ( some chords may be too long )
down=find(xcont>xedge & ycont<-yedge);
up=find(xcont>xedge & ycont>yedge);
cd=polyfit(xcont(down),ycont(down),1);
cu=polyfit(xcont(up),ycont(up),1);
iu1=find(xchord(1,:)>xedge & ychord(1,:)>0 & dum_vangle==-90 );
if ~isempty(iu1)
xchord(1,iu1)=-(b(iu1)-cu(2))./(m(iu1)-cu(1)+eps);
ychord(1,iu1)=m(iu1).*xchord(1,iu1)+b(iu1);
end
iu2=find(xchord(2,:)>xedge & ychord(2,:)>0 & ychord(1,:) & ....
dum_vangle==-90);
if ~isempty(iu2)
xchord(2,iu2)=-(b(iu2)-cu(2))./(m(iu2)-cu(1)+eps);
ychord(2,iu2)=m(iu2).*xchord(2,iu2)+b(iu2);
end
id1=find(xchord(1,:)>xedge & ychord(1,:)<0 & dum_vangle==90);
if ~isempty(id1)
xchord(1,id1)=-(b(id1)-cd(2))./(m(id1)-cd(1)+eps);
ychord(1,id1)=m(id1).*xchord(1,id1)+b(id1);
end
id2=find(xchord(2,:)>xedge & ychord(2,:)<0 & dum_vangle==90);
if ~isempty(id2)
xchord(2,id2)=-(b(id2)-cd(2))./(m(id2)-cd(1)+eps);
ychord(2,id2)=m(id2).*xchord(2,id2)+b(id2);
end
ilow=find(ychord(1,:)<ymin);
ihig=find(ychord(1,:)>ymax);
ilef=find(xchord(1,:)<xmin);
irig=find(xchord(1,:)>xmax);
if ~isempty(ilow)
ychord(1,ilow)=ymin*ones(size(ilow));
xchord(1,ilow)=ymin./m(ilow)-b(ilow)./m(ilow);
end
if ~isempty(ihig)
ychord(1,ihig)=ymax*ones(size(ihig));
xchord(1,ihig)=ymax./m(ihig)-b(ihig)./m(ihig);
end
if ~isempty(ilef)
xchord(1,ilef)=xmin*ones(size(ilef));
ychord(1,ilef)=m(ilef)*xmin+b(ilef);
end
if ~isempty(irig)
xchord(1,irig)=xmax*ones(size(irig));
ychord(1,irig)=m(irig)*xmax+b(irig);
end
ilow=find(ychord(2,:)<ymin);
ihig=find(ychord(2,:)>ymax);
ilef=find(xchord(2,:)<xmin);
irig=find(xchord(2,:)>xmax);
if ~isempty(ilow)
ychord(2,ilow)=ymin*ones(size(ilow));
xchord(2,ilow)=ymin./m(ilow)-b(ilow)./m(ilow);
end
if ~isempty(ihig)
ychord(2,ihig)=ymax*ones(size(ihig));
xchord(2,ihig)=ymax./m(ihig)-b(ihig)./m(ihig);
end
if ~isempty(ilef)
xchord(2,ilef)=xmin*ones(size(ilef));
ychord(2,ilef)=m(ilef)*xmin+b(ilef);
end
if ~isempty(irig)
xchord(2,irig)=xmax*ones(size(irig));
ychord(2,irig)=m(irig)*xmax+b(irig);
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment