diff --git a/TCV/VsxrTCVradius.m b/TCV/VsxrTCVradius.m new file mode 100644 index 0000000000000000000000000000000000000000..38f1ce9f0cd0ec1b071a7d59612bdc9f844bcba9 --- /dev/null +++ b/TCV/VsxrTCVradius.m @@ -0,0 +1,24 @@ +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 diff --git a/TCV/ece_te.m b/TCV/ece_te.m new file mode 100755 index 0000000000000000000000000000000000000000..32f05aa0a82665bd597661807c0bae1a3e3bf9c4 --- /dev/null +++ b/TCV/ece_te.m @@ -0,0 +1,67 @@ +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); + + diff --git a/TCV/get_xtomo_data.m b/TCV/get_xtomo_data.m new file mode 100755 index 0000000000000000000000000000000000000000..3652dd33b3a34c864befe2a31d132a9cd5cbc8b0 --- /dev/null +++ b/TCV/get_xtomo_data.m @@ -0,0 +1,497 @@ +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 + + + + + + diff --git a/TCV/private/angular_fact_1.mat b/TCV/private/angular_fact_1.mat new file mode 100755 index 0000000000000000000000000000000000000000..85fd9a52692c72745b87ad324e2601e6d90fa66d Binary files /dev/null and b/TCV/private/angular_fact_1.mat differ diff --git a/TCV/private/bad_channels.mat b/TCV/private/bad_channels.mat new file mode 100755 index 0000000000000000000000000000000000000000..e950698cfd82fae24b1fcd0aeaf8eea96ef638b1 Binary files /dev/null and b/TCV/private/bad_channels.mat differ diff --git a/TCV/private/ece_d.m b/TCV/private/ece_d.m new file mode 100755 index 0000000000000000000000000000000000000000..17752b9d7456b19dfc7221514fc3f6000d2710f8 --- /dev/null +++ b/TCV/private/ece_d.m @@ -0,0 +1,51 @@ +% 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 ]; + diff --git a/TCV/private/tcv_vesc1.mat b/TCV/private/tcv_vesc1.mat new file mode 100644 index 0000000000000000000000000000000000000000..becc16ea4e9dc15a3afedbacd112ef4f681b7bcb Binary files /dev/null and b/TCV/private/tcv_vesc1.mat differ diff --git a/TCV/xtomo_geometry.m b/TCV/xtomo_geometry.m new file mode 100644 index 0000000000000000000000000000000000000000..a809b9e65cca29d8df8eff407eaac96d7b7ebc4c --- /dev/null +++ b/TCV/xtomo_geometry.m @@ -0,0 +1,402 @@ +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 + + + + +