wave_distribution_yp

 Split the M rays in each wavestructs into N*M rays with P/N power for averaging the effect of plasma fluctuations on ray propagation (C3PO code only). 
 May use distributed processing for reducing the overall computational time

 INPUT
    - dkepath structure
    - dkeparam structure
    - wavestructs structure
    - f0struct structure
    - dt_fluct fluctuation correlation time (is s).

 OUTPUT:
    - waves structure
    - wavestructs structure (it is output since the phase has changed)

 By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)

0001 function [waves,wavestructs,byisel] = wave_distribution_yp(dkepath,dkeparam,wavestructs,f0struct,dt_fluct,byPnonabs)
0002 %
0003 % Split the M rays in each wavestructs into N*M rays with P/N power for averaging the effect of plasma fluctuations on ray propagation (C3PO code only).
0004 % May use distributed processing for reducing the overall computational time
0005 %
0006 % INPUT
0007 %    - dkepath structure
0008 %    - dkeparam structure
0009 %    - wavestructs structure
0010 %    - f0struct structure
0011 %    - dt_fluct fluctuation correlation time (is s).
0012 %
0013 % OUTPUT:
0014 %    - waves structure
0015 %    - wavestructs structure (it is output since the phase has changed)
0016 %
0017 % By Yves Peysson (CEA-DRFC, yves.peysson@cea.fr) and Joan Decker (CEA-DRFC, joan.decker@cea.fr)
0018 %
0019 if nargin < 3,
0020     error('Wrong number of input arguments in wave_distribution_yp.m.');
0021 end
0022 %
0023 if nargin < 6,
0024     byPnonabs = NaN;
0025 end
0026 %
0027 if nargin < 5,
0028     dt_fluct = Inf;
0029 end
0030 %
0031 if nargin < 4,
0032     f0struct = [];
0033 end
0034 %
0035 if isfield(dkeparam,'Nfluct'),
0036     Nfluct = max([1,round(dkeparam.Nfluct)]);
0037 else
0038     Nfluct = 1;
0039 end
0040 %
0041 nw = size(wavestructs,2);
0042 %
0043 % If Nfluct > 1, parallel fluctuation states are considered simultaneously
0044 %
0045 if ~iscell(byPnonabs) && ~isinf(dt_fluct),% first time step, with fluctuations
0046     %
0047     if size(wavestructs,1) > 1,
0048         error('byPnonabs must be provided if wavestruct structure already duplicated.')
0049     end
0050     %
0051     wavestructs = repmat(wavestructs,[Nfluct,1]);% copy wavestructs in Nfluct structures
0052     byPnonabs = cell(Nfluct,nw);
0053     %
0054     for ifluct = 1:Nfluct,
0055         %
0056         for iw = 1:nw,% for each wavestruct
0057             %
0058             % initialize power in each ray by deviding total power
0059             %
0060             if isfield(wavestructs{ifluct,iw}.launch,'bPlhtot'),%LH wave
0061                 %
0062                 wavestructs{ifluct,iw}.launch.bPlhtot_in = wavestructs{ifluct,iw}.launch.bPlhtot/Nfluct;%Total launched LH power is constant (sensitive to core and edge fluctuations)
0063                 byPnonabs{ifluct,iw} = zeros(1,length(wavestructs{ifluct,iw}.launch.bPlhtot));
0064                 %
0065             elseif isfield(wavestructs{ifluct,iw}.launch,'yP_L'),%EC wave
0066                 %
0067                 wavestructs{ifluct,iw}.launch.yP_L_in = wavestructs{ifluct,iw}.launch.yP_L/Nfluct;%Total launched EC power is constant (sensitive to core fluctuations only)
0068                 byPnonabs{ifluct,iw} = zeros(1,length(wavestructs{ifluct,iw}.launch.yP_L));
0069                 %
0070             end
0071             %
0072             % limit tfinal according to dtn_fluct
0073             %
0074             if ~isinf(dt_fluct),
0075                 wavestructs{ifluct,iw}.rayparam.tfinal = dt_fluct*wavestructs{ifluct,iw}.launch.omega_rf;%ray tracing duration cannot exceed the fluctuation correlation time (x 2.3 approximately)
0076             end
0077             %
0078         end
0079     end
0080 end
0081 %
0082 % byPnonabs contains the power not absorbed
0083 %
0084 if iscell(byPnonabs),%calculate fluctuation phases and adjust initial power in each ray
0085     %
0086     if size(wavestructs,1) ~= Nfluct,
0087         error('wavestruct structure must be already duplicated if byPnonabs is provided.')
0088     end
0089     %
0090     for ifluct = 1:Nfluct,
0091         %
0092         if isfield(wavestructs{ifluct,1}.equil_fit,'fluct_fit') && ~isempty(fieldnames(wavestructs{ifluct,1}.equil_fit.fluct_fit)),
0093             %
0094             fluct_fit = fluctphase_yp(wavestructs{ifluct,1}.equil_fit.fluct_fit);%Calculations with new random phase for each structure (equilibrium fluctuations). Spectrum fluctuations are done in main_rayinit_launch_jd.m
0095             %
0096             for iw = 1:nw,
0097                 %
0098                 wavestructs{ifluct,iw}.equil_fit.fluct_fit = fluct_fit;%Replicate the same fluctuation structure to all wavestructs (equilibrium fluctuations)
0099                 %
0100             end
0101         end
0102         %
0103         for iw = 1:nw,
0104             if isfield(wavestructs{ifluct,iw}.launch,'bPlhtot'),%LH wave
0105                 %
0106                 wavestructs{ifluct,iw}.launch.bPlhtot = wavestructs{ifluct,iw}.launch.bPlhtot_in + byPnonabs{ifluct,iw};%Total launched LH power is constant
0107                 %
0108             elseif isfield(wavestructs{iw}.launch,'yP_L'),%EC wave
0109                 %
0110                 wavestructs{ifluct,iw}.launch.yP_L = wavestructs{ifluct,iw}.launch.yP_L_in + byPnonabs{ifluct,iw};%Total launched LH power is constant
0111                 %
0112             end          
0113         end 
0114     end          
0115 end    
0116 %
0117 dkecluster = clustermode_luke(dkeparam.clustermode,'wave_solver_yp',dkepath);%MatLab distributed computing environment
0118 %
0119 if dkecluster.scheduler.mode <= 2,
0120     [flag,waves] = mdce_run(@wave_solver_yp,{dkepath,'',f0struct},2,wavestructs(:),dkecluster);
0121 else% remote computing
0122     [flag,waves] = mdce_run(@wave_solver_yp,{'','',f0struct},2,wavestructs(:),dkecluster);
0123 end    
0124 %
0125 % determines which rays actually enter the plasma and report the ray number in zP_0_2piRp_mod
0126 %
0127 byisel = cell(Nfluct,nw);
0128 iray = 0;
0129 for iw = 1:nw,
0130     for ifluct = 1:Nfluct,
0131         wave = waves{ifluct+Nfluct*(iw-1)};
0132         %
0133         if isfield(wavestructs{ifluct,iw}.launch,'bPlhtot'),%LH wave
0134             ny = length(wavestructs{ifluct,iw}.launch.bPlhtot);%number of lobes
0135         elseif isfield(wavestructs{iw}.launch,'yP_L'),%EC wave
0136             ny = length(wavestructs{ifluct,iw}.launch.yP_L);%number of rays
0137         end          
0138         %
0139         byisel{ifluct,iw} = cell(1,ny);
0140         %
0141         for iy = 1:ny,
0142             byisel{ifluct,iw}{iy} = [];
0143         end
0144         %
0145         % byisel contains the ray number belonging to the same lobe with possibly different launching positions (LH) or
0146         % to the same beam with different beamlets (EC)
0147         %
0148         for ib = 1:length(wave.rays),
0149             iray = iray+1;
0150             %
0151             iy = wave.rays{ib}.rayinits{1}.yib;
0152             %
0153             byisel{ifluct,iw}{iy} = [byisel{ifluct,iw}{iy},iray];
0154         end
0155     end
0156 end
0157 %