make_wavelaunchs_metis2luke_ec

 This function returns EC wave launching parameters

 by J. Decker and J.-F. Artaud, 02/02/2011

0001 function launchs = make_wavelaunchs_metis2luke_ec(option,ecparam,equil)
0002 %
0003 % This function returns EC wave launching parameters
0004 %
0005 % by J. Decker and J.-F. Artaud, 02/02/2011
0006 %
0007 if ecparam.pec < ecparam.pecmin,% not enough power to account for wave
0008     launchs = cell(1,0);
0009     return
0010 end
0011 %
0012 [qe,me,void,void,void,void,void,mc2] = pc_dke_yp;
0013 %
0014 offset = 0.01;
0015 %
0016 if ecparam.mode == 1,
0017     %
0018     % C3PO/LUKE case with one launcher and horizontal propagation. The
0019     % toroidal angle is chosen to match the prescribed poloidal location
0020     %
0021     launch.id = 'EC';
0022     launch.type = 'EC';
0023     %
0024     launch.dNpar0 = ecparam.dnpar0;%ray spectral width for Fokker-Planck calculations
0025     launch.mmode = ecparam.mmode;% (-1) for O mode (1) for X mode
0026     %
0027     launch.ns = ecparam.ns;
0028     launch.method = ecparam.method;
0029     %
0030     launch.yP_L = ecparam.pec;% wave power (W)
0031     %
0032     nr = length(equil.psi_apRp);
0033     %
0034     px = linspace(0,1,nr);% metis radial grid
0035     %
0036     xec = ecparam.xec;% radial position at deposition location
0037     polec = ecparam.polec*pi/180;% poloidal angle at deposition location (rad)
0038     omega_rf = ecparam.freqec*2*pi*1e9;%wave angular frequency (rad/s)
0039     alpha_L = ecparam.alpha_L*pi/180;% toroidal launching angle (rad)
0040     %
0041     pabs = ecparam.pabs;%momentum at absorption (for Doppler shift)
0042     %
0043     if isnan(polec) && isnan(omega_rf),
0044         error('One parameter among these two must be specified : poloidal location and EC frequency')
0045     end
0046     if ~isnan(polec) && ~isnan(alpha_L) && omega_rf > 0,
0047         disp('Warning : poloidal location, EC frequency, and toroidal launching angle cannot be specified independently. polodal location discarded')
0048         polec = NaN;
0049     end
0050     %
0051     R_L = 1.1*max(max(equil.Rp + equil.ptx));
0052     beta_L = pi/2;%horizontal propagation
0053     phi_L = ecparam.phi0;
0054     %
0055     yx_L = R_L - equil.Rp;
0056     xz_L = phi_L*equil.Rp;
0057     %
0058     % straight line to reach virtual launching point (on a vertical line
0059     % located at R = 1.1 Rmax where Rmax is the maximum R in the plasma)
0060     %
0061     s = linspace(0,2*equil.Rp,launch.ns);
0062     %
0063     Te_abs = interp1(px,equil.pTe,ecparam.xec,ecparam.method);
0064     betath = sqrt(Te_abs/mc2);
0065     %
0066     if ~isnan(polec),%Doppler shift used to determine frequency or launching angle
0067         %
0068         R_abs = equil.Rp + interp2(px,equil.theta,equil.ptx.',xec,polec,ecparam.method);
0069         Z_abs = equil.Zp + interp2(px,equil.theta,equil.pty.',xec,polec,ecparam.method);
0070         Bx_abs = interp2(px,equil.theta,equil.ptBx.',xec,polec,ecparam.method);
0071         By_abs = interp2(px,equil.theta,equil.ptBy.',xec,polec,ecparam.method);
0072         BPHI_abs = interp2(px,equil.theta,equil.ptBPHI.',xec,polec,ecparam.method);
0073         %
0074         B_abs = sqrt(Bx_abs^2 + By_abs^2 + BPHI_abs^2);
0075         %
0076         Z_L = Z_abs;%horizontal propagation
0077         %
0078         if R_abs < equil.Rp && abs(Z_L - equil.Zp)/max(equil.pty(:)) < offset,
0079             disp(['Warning : Propagation right across the center is problematic in C3PO. Small relative vertical shift of ',num2str(offset),' enforced.'])
0080             Z_L = Z_L + offset*max(equil.pty(:));
0081         end
0082         %
0083         yx_abs = R_abs - equil.Rp;
0084         yy_abs = Z_abs - equil.Zp;
0085         yy_L = Z_L - equil.Zp;
0086         %
0087         if omega_rf > 0,%alpha from Doppler shift
0088             %
0089             B1_res = omega_rf*me/qe; 
0090             n_ec = fix(B1_res/B_abs);%harmonic number
0091             deltalnR = B1_res/(n_ec*B_abs) - 1;%relative B-field gap
0092             %
0093             if deltalnR > 1/(2*n_ec),
0094                 disp('Warning : prescribed absorption appears to be behind nominal resonance, alpha = 0 enforced')
0095                 sinalpha_abs = 0;
0096             else
0097                 sinalpha_abs = deltalnR/pabs/betath;%approximate Doppler shift estimation, based on Decker PoP 06 with pnp = 2
0098                 if sinalpha_abs > 1.0,
0099                     disp('Warning : prescribed Doppler shift appears to be too large, local alpha limited to pi/2')
0100                     sinalpha_abs = 1.0;
0101                 end
0102             end
0103             %
0104             alpha_abs = option.signe*option.sens*asin(sinalpha_abs);
0105             %
0106             [void,void,syR,void,void,syalpha] = raypath_prop_jd(equil.Rp,yx_abs,yy_abs,xz_L,alpha_abs,beta_L,s);
0107             %
0108             alpha_L = pi + interp1(syR,syalpha,R_L,'linear');
0109             %
0110         else%omega from Doppler shift
0111             %
0112             n_ec = -omega_rf/(2*pi*1e9);
0113             %
0114             if ~isnan(alpha_L),
0115                 %
0116                 [void,void,syR,void,void,syalpha] = raypath_prop_jd(equil.Rp,yx_L,yy_L,xz_L,alpha_L,beta_L,s);
0117                 %
0118                 alpha_abs = interp1(syR,syalpha,R_abs,'linear');
0119                 %
0120                 if isnan(alpha_abs),
0121                     disp('Warning : Doppler shift could not be determined, enforced to 0')
0122                     alpha_abs = 0;
0123                 end
0124                 %
0125             else % alpha_L determined from local alpha
0126                 %
0127                 if option.sens ~= 0,% ECCD
0128                     %
0129                     alpha_abs = option.signe*option.sens*30*pi/180;%approximate optimal angle
0130                     %
0131                 else                % ECRH
0132                     %
0133                     alpha_abs = 0;
0134                     %
0135                 end                    
0136                 %
0137                 [void,void,syR,void,void,syalpha] = raypath_prop_jd(equil.Rp,yx_abs,yy_abs,xz_L,alpha_abs,beta_L,s);   
0138                 %
0139                 alpha_L = pi + interp1(syR,syalpha,R_L,'linear');                    
0140                 %
0141             end
0142             %
0143             sinalpha_abs = abs(sin(alpha_abs));
0144             %
0145             deltalnR = sinalpha_abs*pabs*betath;%approximate Doppler shift estimation, based on Decker PoP 06 with pnp = 2
0146             %
0147             B_res = B_abs*(1 + deltalnR);%relative B-field gap
0148             %
0149             omega_rf = n_ec*B_res*qe/me;
0150             %
0151         end
0152         %
0153     else% deposition location from wave characteristics
0154         %
0155         tR_abs = equil.Rp + interp1(px,equil.ptx,xec,ecparam.method);
0156         tZ_abs = equil.Zp + interp1(px,equil.pty,xec,ecparam.method);
0157         tBx_abs = interp1(px,equil.ptBx,xec,ecparam.method);
0158         tBy_abs = interp1(px,equil.ptBx,xec,ecparam.method);
0159         tBPHI_abs = interp1(px,equil.ptBPHI,xec,ecparam.method);
0160         %
0161         [void,void,syR,void,void,syalpha] = raypath_prop_jd(equil.Rp,yx_L,0,xz_L,alpha_L,beta_L,s);
0162         %
0163         talpha_abs = interp1(syR,syalpha,tR_abs,'linear');
0164         %
0165         tsinalpha = abs(sin(talpha_abs));
0166         %
0167         tdeltalnR = tsinalpha*pabs*betath;%approximate Doppler shift estimation, based on Decker PoP 06 with pnp = 2
0168         %
0169         tB_abs = sqrt(tBx_abs.^2 + tBy_abs.^2 + tBPHI_abs.^2);
0170         %
0171         B1_res = omega_rf*me/qe; 
0172         n_ec = max(fix(B1_res./tB_abs));%harmonic number
0173         %
0174         tB1_res = n_ec.*tB_abs.*(1 + tdeltalnR);%relative B-field gap
0175         %
0176         tmask = 1:(length(tB1_res) + 1)/2;%upper midplane
0177         %
0178         Z_L = interp1(tB1_res(tmask),tZ_abs(tmask),B1_res);
0179         %
0180         if isnan(Z_L),
0181             disp('Warning : optimum Z_L not found, Z_L = 0 enforced')
0182             Z_L = 0;
0183         end            
0184         %
0185     end
0186     %
0187     launch.omega_rf = omega_rf;
0188     %
0189     if alpha_L > pi,
0190         alpha_L = alpha_L - 2*pi;
0191     end
0192     %
0193     if option.signe*option.sens*alpha_L > 0,
0194         disp('Warning : calculated launching angles not compatible with co/counter current prescrion, hence modified')
0195         alpha_L = -alpha_L;
0196     end
0197     %
0198     launch.yR_L = R_L;
0199     launch.yZ_L = Z_L;
0200     launch.yalpha_L = alpha_L;
0201     launch.ybeta_L = beta_L;
0202     launch.yphi_L = phi_L;
0203     %
0204     launchs{1} = launch;
0205     %
0206 elseif ecparam.mode == 2,
0207     %
0208     % C3PO/LUKE case with detailed launcher structure
0209     %
0210     error('ecparam.mode = 2 not yet implemented')
0211     %
0212 end
0213 %