Description of reflect

0001 function [P_ref, K_ref, E_ref] = reflect_tcv(shotnum, shotime , P_end, K_end, E_end )
0002 % ray freq. X2-mode in TCV 82.7 GHz
0003 %
0004 %
0005 %   reflect_tcv.m
0006 %
0007 % inputs: (from LUKE/C3PO LCFS)
0008 % - shot and time numbersS
0009 % - X exit ray position from LCFS (p) in [m]
0010 % - K wave vector in X in [m-1], gives also the direction of the ray propagation
0011 % - E field in X
0012 %
0013 % X,K,E are in the rectangurlar reference frame of TCV (used by TORAY)
0014 %
0015 %       - x is the major radius direction: 0 at the axis of the machine; positive
0016 %         outward, from Left to Right in a TCV top view)
0017 %       - y is the major radius direction: 0 at the axis of the machine; positive
0018 %         upward, from the axis in a TCV top view)
0019 %       - z is the vertical direction: 0 at the axis of the machine; positive
0020 %         coming out from the axis in a TCV top view)
0021 %
0022 % to interface with LUKE/C3PO reference (rho,theta,phi) we need to tranform the
0023 % coodinates
0024 %re
0025 % outputs:
0026 % - WaveX ray component in X-MODE (XoutX,KoutX,EoutX)
0027 %   NOTE: the Kout is normalized to k_in ( kk = -k_in/norm(k_in) );
0028 %   thus, norm(Kout) = 1
0029 %
0030 % - WaveO ray component in O-MODE (XoutO,KoutO,EoutO)
0031 % - P contains the power fraction in X and O-MODE in the new incoming LCFS position
0032 %
0033 % the main body of this function exploit the reflect_toray_rays.m routine of
0034 %  /home/ffelici/matlab/
0035 
0036 
0037 addpath /home/ffelici/matlab/ECPOL/pol_code_develop/
0038 addpath /home/ffelici/matlab/ECPOL/pol_code_develop/ECPOL/
0039 addpath /home/gnesin/
0040 
0041 %
0042 R_end_XY = sqrt(P_end(1).^2+ P_end(2).^2);
0043 
0044 
0045 
0046 
0047 
0048 X=P_end; K = K_end; E = E_end;       
0049 
0050 disp(' ');
0051 disp('*** running reflect_tcv ***')
0052 disp(' ');
0053 
0054 shot = str2num(shotnum);
0055 time = str2num(shotime);
0056  %shot = 25474; time = 0.9;
0057      
0058 do_polarization = 1; % do not worry about polarization for now, just flip k vector
0059 doplot = 1;
0060 plot_reflection = 0;
0061 do_top_plot =0; % =1 enable 3D plot to check the ray reflection
0062 
0063 % data at plasma exit (from LUKE/C3PO) the calculation is performed ray by ray
0064  % the Xposition is in [m]
0065   
0066  %iend(iray) = find(~isnan(raydata.s_out(:,iray)),1,'last'); %index upon which ray ends
0067  xend = X(1); %(:,iray) = raydata.x_out(iend(iray),iray)/100;
0068  yend = X(2); %(:,iray) = raydata.y_out(iend(iray),iray)/100;
0069  zend = X(3); %(:,iray) = raydata.z_out(iend(iray),iray)/100;
0070  %pend(:,iray) = raydata.p_out(iend(iray),iray)/100;
0071  kxend= K(1); %(:,iray) = raydata.kvx(iend(iray),iray);
0072  kyend= K(2); %(:,iray) = raydata.kvy(iend(iray),iray);
0073  kzend= K(3); %(:,iray) = raydata.kvz(iend(iray),iray);
0074  
0075  Exend= E(1); %(:,iray) = [raydata.ex_r(iend(iray),iray) + i*raydata.ex_i(iend(iray),iray)];
0076  Eyend= E(2); %(:,iray) = [raydata.ey_r(iend(iray),iray) + i*raydata.ey_i(iend(iray),iray)];
0077  Ezend= E(3); %(:,iray) = [raydata.ez_r(iend(iray),iray) + i*raydata.ez_i(iend(iray),iray)];
0078  
0079 % position, k vectors and E field at plasma exit
0080 %pend = [1.04;-0.45;0];% just for test
0081 pend = [xend;yend;zend]; % position
0082 %kend = [1.0834;-0.3321;-0.0006].*(1000); % just for test
0083 kend = [kxend;kyend;kzend]; % k vector
0084 Eend = [Exend;Eyend;Ezend]; % E vector
0085 
0086 
0087 
0088 % check the end point is in the TCV vacuum chamber
0089 
0090 
0091 %----------------------------------------------------------
0092 disp(' '), disp(' ray parameters (exit from the LCFS): (RED point)')
0093 disp(['P_exit = ', num2str(pend')])
0094 disp(['k_exit = ', num2str(kend')])
0095 disp(['E_exit = ', num2str(Eend')]), 
0096 
0097 
0098 
0099 
0100 me = 9.10938188e-31; qe = 1.60217646e-19;
0101 %O_c = qe*sqrt(sum(B.^2,1)) / me / (2*pi);
0102 freq = 82.7e9;
0103 cs = 299792458; % speed of light [m/s]
0104 
0105 kkk = (freq*2*pi)/cs; % test the k norm
0106 disp(['NORM(K_test) = freq/c: ', num2str(kkk), ' [1/m]']),disp(' ')
0107 
0108 %------------------- intersect with the TCV wall -----------
0109 mdsopen(shot)
0110 rt=mdsdata('static("R_T")'); % in cm!
0111 zt=mdsdata('static("Z_T")');
0112 rr=rt;
0113 zz=zt;  
0114 
0115 
0116 % allocate
0117 inointers = [];
0118 l_out = nan; %*ones(raydata.nbr,1);
0119 xw = nan; %*ones(raydata.nbr,1);
0120 yw = nan; %*ones(raydata.nbr,1);
0121 zw = nan; %*ones(raydata.nbr,1);
0122 
0123 
0124 % here use inters3d to find the distance from the exit point to the TCV wall
0125  [l2pl,dd,cc]=inters3d_tcv(xend,yend,zend,...
0126                        kxend,kyend,kzend,rr',zz');
0127 
0128  %if  length(l2pl)<=2 & (R_end_XY > 1.13596 | R_end_XY < 0.624 |P_end(3) < -0.75 | P_end(3) > 0.75)
0129  %l_out = max(abs(l2pl)); % to consider P_out on the wall
0130  %disp(' WARNING : l_out is chosen as a "max"')
0131  %else
0132  %l_out = min(abs(l2pl)); % distance of closest wall intersection
0133  l_out =min(l2pl(find(l2pl>100.*eps)));
0134  
0135  if sum(sign(l2pl)) == - length(l2pl)
0136  disp('l2pl has only negative elements')
0137  l_out = 0,min(abs(l2pl(find(abs(l2pl)>100.*eps))));
0138  end
0139  
0140  
0141  %end;
0142  
0143  if ~isempty(l_out)%(Pw),            % if 0 then no intersection points
0144                                         % wall intersection points
0145                                         xw= xend+kxend*l_out;%Pw(1); %xend+kxend*l_out;
0146                                         yw=yend+kyend*l_out;%Pw(2); %yend+kyend*l_out;
0147                                         zw=zend+kzend*l_out;%Pw(3); %zend+kzend*l_out;
0148  else
0149   disp('no intersection found with the wall');
0150   disp('intersection point is replaced with the edge point')
0151   inointers = [inointers 1]; % indices of rays with no wall intersection
0152  end
0153  
0154   if ~isempty(inointers)
0155    warning(['****** ', int2str(numel(inointers)), ' ray(s) did not find an intersection, eliminating ray(s)']) 
0156   end 
0157 
0158 rw=sqrt(xw.^2+yw.^2);        % major-radial location of intersection point
0159 inotinnerwall = find((rw-min(rr))>1e-6); % rays not hitting inner wall
0160 rw1 = (xw.^2+yw.^2+zw.^2);
0161 % ----------------------------------------------------------------------
0162 %    REFLECTION CASES; TCV CROSS SECTION IS NOT EXACLTLY A RECTANGLE!!!
0163 % ----------------------------------------------------------------------
0164 R_tcv_int = 0.624 + 0.0001; R_tcv_out = 1.136 - eps;
0165 direction = 0;
0166 sign_w =1; % it is ok for outer reflection (also top and bottom)
0167 
0168  if abs(xw)  <= R_tcv_int & abs(yw)  <= R_tcv_int
0169     sign_w = -1;
0170     disp('Inward ray: inner wall reflection.')
0171  
0172     
0173  end; 
0174  
0175  if abs(rw)  >= R_tcv_int & abs(zw)  <= 0.55
0176     sign_w=1  ;  
0177     disp(' Outward ray: outer wall reflection')
0178  end
0179  
0180  if abs(rw)  >= R_tcv_int & (zw)  >= (max(zz)-0.0001)
0181     sign_w=1  ;  
0182     direction = 1;
0183     disp(' TOP wall reflection')
0184  end
0185  
0186  if abs(rw1) > (0.97.^2  + 0.75.^2  ) & abs(rw1) < (1.135.^2  + 0.55.^2  ) & zw>0.55
0187     sign_w=1  ;  
0188     direction =  2;
0189     th_2 = atan((0.75-0.55)./(1.135-0.97));  % angle for the normal vector in the X-Z plane
0190     phi_2 = atan((yw)./(xw));
0191     disp(' TOP(out) wall reflection')
0192  end
0193   
0194  if abs(rw1) > (0.62457.^2  + 0.704.^2  ) & abs(rw1) < (0.67164.^2  + 0.75.^2  ) & zw>0.62457
0195     sign_w=1  ;  
0196     direction =  3;
0197     th_3 = atan((0.75- 0.704)./(0.67164-0.62457));% is 45 degrees
0198     phi_3 = atan((yw)./(xw));
0199     disp(' TOP(in) wall reflection')
0200  end
0201  
0202  if abs(rw)  >= R_tcv_int & (zw)  <= (min(zz)+0.0001)
0203     sign_w=1  ;  
0204     direction =  -1;
0205     disp(' BOTTOM wall reflection')
0206  end
0207 
0208 
0209 if abs(rw1) > (0.97.^2  + 0.75.^2  ) & abs(rw1) < (1.135.^2  + 0.55.^2  ) & zw<-0.55
0210     sign_w=1  ;  
0211     direction =  -2;
0212     th_2 = atan((0.75-0.55)./(1.135-0.97));
0213     phi_2 = atan((yw)./(xw));
0214     disp(' BOTTOM(out) wall reflection')
0215  end
0216 
0217  
0218  if abs(rw1) > (0.62457.^2  + 0.704.^2  ) & abs(rw1) < (0.67164.^2  + 0.75.^2  ) & zw<0.62457
0219     sign_w=1  ;  
0220     direction =  -3;
0221     th_3 = atan((0.75- 0.704)./(0.67164-0.62457)); % is 45 degrees
0222     phi_3 = atan((yw)./(xw));
0223     disp(' BOTTOM(in) wall reflection')
0224  end
0225 
0226 % --------------------------------------------------
0227 
0228 if ~isempty(inotinnerwall)
0229  warning(['****** ', int2str(numel(inotinnerwall)), ' ray(s) did not intersect inner wall, eliminating ray(s)']) 
0230 end
0231 
0232 %iok = setdiff(i_notelim,union(inointers,inotinnerwall)); % iok takes the index of the remaining
0233 % beams over which the calculation will go on.. i = index ok = that are OK for the prosecution
0234 % of the current analysis
0235 
0236 %----------------------------------
0237  iok = 1; % for a ray by ray case
0238 %----------------------------------
0239 
0240 if isempty(iok); 
0241  warning('no rays could be reflected')
0242  raygeom_files = []; torinput_files = [],suffixes=[]; % return empty outputs
0243  return
0244 end
0245 
0246 %----------------------- inward Wall normals at intersection points--
0247 zw = zw; % z of wall
0248 %sign_w= 1;
0249 nw = sign_w.*[xw'; yw'; zeros(size(xw))'] ./ (ones(3,1)*rw'); % inner wall normal in the (X-Y)plane
0250 pw = [xw';yw';zw']; % position of wall
0251 
0252 disp('the Black point indicates the'),disp([' wall interctepion at p_wall = ' num2str(pw')])
0253 
0254  % Use ECPOL tools to handle the reflection including the correct
0255 % polarization
0256 if ~exist('mirror','file'); error('ECPOL tools are not in your path'); end
0257 
0258 k_in = kend(:,iok); % "in" means INTO the mirror (wall)
0259 E_in = Eend(:,iok);
0260 
0261 % if we don't worry about the polarization, any E field will do (ignored
0262 % anyway)
0263 
0264 %---------------------------------------------------
0265 if do_polarization
0266     % pre-load some stuff we'll need
0267     %load(fname_torinputs)
0268     %time = torinputs.time;
0269     %shot = torinputs.shot;
0270     
0271     psi = psitbxtcv(shot,time);  % get psi from the tree. This is psi with the correct signs and magnitude at the center.
0272     
0273     mdsopen(shot);
0274     rbou = mdsdata('\results::r_contour[,$1]',time); rbou = rbou(~isnan(rbou));
0275     zbou = mdsdata('\results::z_contour[,$1]',time); zbou = zbou(~isnan(zbou));
0276     
0277     % (rbou,zbou)are the coordinates of the last close flux surface from the TCV
0278     %  magnetic reconstruction
0279     
0280     RBphi = mdsdata('\magnetics::RBPHI[$1]',time); % get R*Bphi (is a scalar)
0281     
0282     mdsclose;
0283 end
0284 %--------------------------------------------------------------------------
0285 %------------------
0286 %-----------------------------------------------
0287 %  ANALISI AT THE EXIT OF THE LCFS P = P_end
0288 % -----------------------------------------------
0289  if ~do_polarization
0290     Nk = null(k_in'); E_in = Nk(:,1); % E is perpendicular to k
0291     wavein = wave(k_in,E_in,raygeomdata(iok).freq*1e9);
0292   else    
0293     % determine local B field - copied from wiki psitbx page
0294     x=pend(1); y=pend(2); z=pend(3);
0295     pgrid = psitbxgrid('orthogonal','points',{x,y,z}); % define point at which we want the B
0296     dpsidz = psitbxf2f(psi,pgrid,[0 1 0]); % spatial derivatives if psi
0297     dpsidr = psitbxf2f(psi,pgrid,[1 0 0]);
0298     R = sqrt(x.^2+y.^2); % R of the point
0299     Br = -1./R/(2*pi).*dpsidz.x; % Radial component of B
0300     Bz =  1./R/(2*pi).*dpsidr.x; % Z component of B
0301     Bphi = RBphi/R; % Toroidal component of B (positive anticlockwise looking from the top)
0302     % NB this way of computing Bphi is only valid outside the plasma!
0303     % now get the xyz components of the B field
0304     RR = [[x;y],[-y;x]]./norm([x;y]); % transformation matrix to get to xyz from cylindrical coord
0305     Bxy = RR*[Br;Bphi];
0306     
0307     %--------------------------------------------------------
0308     %---------------- B field at the Plasma exit form the LCFS:(P_end)
0309     
0310     Bxyz= [Bxy(1);Bxy(2);Bz];
0311     % all of the above was to get the local B field
0312     
0313     %---------------------------------------------------------
0314     
0315     kk = -k_in/norm(k_in); % normalized k vector back INTO the plasma (- sign)
0316     %kk is the same as (-1).*vk_out
0317     
0318     freq_thisray = 82.7*1e9;
0319     %freq_thisray = raygeomdata(iok).freq*1e9; % frequency of this ray [in Hz]
0320     [Xwave,Owave] = ox_coupling(kk,Bxyz,freq_thisray);
0321     % returns wave corresponding to O and X mode coming OUT of the plasma at LCFS?
0322     % flip k again
0323     Xwave.k = -Xwave.k; % k vector at P_end is the same as vk_out
0324     Owave.k = -Owave.k;  
0325     
0326     Xwave.Position = pend ; % 100*pend;  % check the units
0327     Owave.Position = pend ; % 100*pend;  % here we need [m]
0328   
0329   %--------------------- COMMENTS ------------------------
0330   % the O/X wave structure conteins the k & E vectors of the wave leaving the LCFS
0331   % and directed towards the TCV wall for a reflection
0332   % -------------------------------------------------------------------------
0333   
0334  end 
0335  %%%%%%%%%%%%%%%%%%%%%%%%%% the part here above can be deleted
0336  %%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%
0337  
0338 mywave = wave(k_in,E_in,freq_thisray);
0339 mywave.Position = pw;
0340 
0341 mirr = mirror('flat',nw(:,iok)); % mirror with normal corresponding to wall normal
0342  %mirr = mirror('flat',n_Pw); % mirror with normal corresponding to wall normal
0343 
0344  % ------------------------------------------------------------------
0345  % differents wall orientations requires adapted normal vector (nw)
0346  
0347 if direction == 1;
0348 nw=[0;0;1]; % pw is on the TCV TOP
0349 mirr = mirror('flat',nw); 
0350 end;
0351 
0352 if direction == 2;
0353 nw=[sin(th_2);sin(phi_2);cos(th_2)]; % pw is on the TCV TOP-out
0354 
0355   if xw < 0
0356       nw= [-sin(th_2);-sin(phi_2);cos(th_2)]; %  is ok for X<0, Y<0
0357   end
0358 
0359 mirr = mirror('flat',nw); 
0360 end;
0361 
0362 if direction == 3;
0363 nw=[-cos(pi/4);-sin(phi_3);sin(pi/4)];%/norm([-1;0;1]); % pw is on the TCV TOP-in
0364  if xw<0
0365      nw=[cos(pi/4);sin(phi_3);sin(pi/4)];%
0366   end
0367 
0368 mirr = mirror('flat',nw); 
0369 end;
0370 
0371 if direction == -1
0372 nw= [0;0;-1]; % pw is on the TCV BOTTOM
0373 mirr = mirror('flat',nw); 
0374 end
0375 
0376 if direction == -2
0377 nw= [cos(th_2);sin(phi_2);-sin(th_2)]; % pw is on the TCV BOTTOM-out
0378  if xw < 0
0379       nw= [-cos(th_2);-sin(phi_2);-sin(th_2)]; %
0380   end;
0381 mirr = mirror('flat',nw); 
0382 end
0383 
0384 if direction == -3;
0385 nw=[-cos(pi/4);-sin(phi_3);-sin(pi/4)];%/norm([-1;0;-1]); % pw is on the TCV BOTTOM-in
0386  if xw<0 
0387      nw=[cos(pi/4);sin(phi_3);-sin(pi/4)];
0388   end
0389 
0390 mirr = mirror('flat',nw); 
0391 end;
0392  % -------------------------------------------------------------
0393  if do_polarization
0394    % Reflect waves off wall
0395    % the wave out structure contains the k,E,freq of the wave reflected at the wall
0396    waveoutX = reflect(Xwave,mirr); %
0397    waveoutO = reflect(Owave,mirr);
0398    
0399    waveout = reflect(mywave,mirr);
0400    
0401    % ----------------- COMMENTS ---------------------
0402    % The waveout structure contains the k reflected at the wall;
0403    % -------------------------------------------------
0404    
0405  else
0406     waveoutX = reflect(wavein,mirr);
0407     coupling_matrix = [];
0408  end
0409 
0410 
0411  
0412 Pw=[xw;yw;zw]; vk_out = K_end./norm(K_end); vk_ref2=waveoutX.k; n_Pw = nw;
0413 
0414   if do_top_plot ==1  % plot the ray interception with TCV wall
0415      reflect_plot; % in /home/gnesin/
0416   end;
0417 
0418 % ------------------ Prepare OUTPUTS --------------------------------
0419 P_ref =  Pw; %[xw yw zw]';%= waverefX.Position;  % point of reflection at the TCV wall
0420 K_ref = (waveout.k); % % relected wave number
0421 E_ref = waveout.E; 
0422 
0423 %save /home/gnesin/HXR/Mission_865_X2X3/shot_25474/reflect_ray.mat P_ref K_ref EX_ref EO_ref OX_ref
0424 
0425 disp(' '); disp(' OUTPUTS: K_ref; E_ref; P_ref; ')
0426    disp(['P_ref = ',num2str(P_ref')])
0427    disp(['K_ref = ',num2str(K_ref'),'; norm(K_ref) = 1']) %
0428    disp(['E_ref = ',num2str(E_ref'),' ' ])
0429
reflect_tcv

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
