proc_aloha_spec_jd

 function for coupling ALOHA outputs and C3PO inputs. 
   parametrization can be performed interactively

 INPUTS : 
       - zN_in : initial ntor grid [1,nz]
       - zdN_in : initial ntor grid steps [1,nz] or [1,1] or NaN 
       - zdP_in : initial spectrum dP/dntor [1,nz]
       - Pin_in : corrected initial total power (NaN to read Pin from aloha, 1 for normalized output)
           imaginary part : see param.dir below (used only if param is a scalar)
       - param_in : parametrization of aloha spectrum treatment.
           Note : not all options are available through cases 1 & 2.

         1) Case where param is a structure : fields ->   
           o nlobes : number of lobes selected

           o dir : sort peaks according to : 
                    (0) P     + preserve repartition and total power 
                    (1) P/N^2 + preserve directivity and repartition 
                    (2) P/N^2 + preserve directivity and total power (treating positive and negative spectra as blocks. Warning : some transfer of P/N^2 across N=0 necessary in this case)
                    (3) P/N^2 + preserve repartition and total power 

           o mode : (0) select a number of Gaussian peaks left unchanged (WARNING : total power less than Pin)
                    (1) assemble peaks fixing Pmax and N, allowing dN to change
                    (2) assemble peaks allowing Pmax, N and dN to change

         2) Case where param is empty : to be determined interactively

         3) Case where param is a scalar : param = x + i*y. (default : 2i)
                 The lobes are selected in the order of decreasing power.
                 
           o case x = 0 : lobes calculated by keeping N and Pmax of
               each initial selected lobe. dN increased to account for
               power of neighboring lobes. 
                   - If 0 < y < 1, y is mimimum respective fraction of max power
                   - If y > 1 (integer), y is number of lobes selected
                   - If y < -1 (integer), the power cannot be transfered across N=0

           o case x > 0 : lobes calculated by varying N, Pmax and dN to account for
               neighboring lobes (weighted averages). 
                   - If x < 1, x is mimimum respective fraction of max power
                   - If x > 1 (integer), x is number of lobes selected
                   - y (integer >= 0) is the number of completely protected lobes

           o case x < 0 : identical to x = 0, but at the end only |x| (integer) 
               lobes are truly returned. The lobes are picked alternatively in each direction
             WARNING : in this case the total power is less than Pin
                   
       - opt_optim : use of function lsqnonlin for gaussian lobe fit (default : 0)
       - sdisplay : display structure, with the following field
           o inter : interactive mode for selecting parameters
           o mode : display (0) none (1) text only (2) text and figures (default : 0)
           o style (optional) : style properties (font...) for interactive windows and figures
           o opt_gui (optional) : use gui for interactive windows, use as handle if opt_gui is handle
           o axs (optional) : use axs to plot spectrum (iluke only)
           

 OUTPUTS : 
       - yN : values for output lobes
       - ydN : respective gaussian width
       - yP : respective power
       - param_out : processing parameter used to generate spectrum
       - Pin : initial power used to generate spectrum
       - zyfac : distribution matrix

 by J. Decker and Y. Peysson

0001 function [yN,ydN,yP,param_out,Pin,zyfac] = proc_aloha_spec_jd(zN_in,zdN_in,zdP_in,Pin_in,param_in,opt_optim,sdisplay)
0002     %
0003     % function for coupling ALOHA outputs and C3PO inputs.
0004     %   parametrization can be performed interactively
0005     %
0006     % INPUTS :
0007     %       - zN_in : initial ntor grid [1,nz]
0008     %       - zdN_in : initial ntor grid steps [1,nz] or [1,1] or NaN
0009     %       - zdP_in : initial spectrum dP/dntor [1,nz]
0010     %       - Pin_in : corrected initial total power (NaN to read Pin from aloha, 1 for normalized output)
0011     %           imaginary part : see param.dir below (used only if param is a scalar)
0012     %       - param_in : parametrization of aloha spectrum treatment.
0013     %           Note : not all options are available through cases 1 & 2.
0014     %
0015     %         1) Case where param is a structure : fields ->
0016     %           o nlobes : number of lobes selected
0017     %
0018     %           o dir : sort peaks according to :
0019     %                    (0) P     + preserve repartition and total power
0020     %                    (1) P/N^2 + preserve directivity and repartition
0021     %                    (2) P/N^2 + preserve directivity and total power (treating positive and negative spectra as blocks. Warning : some transfer of P/N^2 across N=0 necessary in this case)
0022     %                    (3) P/N^2 + preserve repartition and total power
0023     %
0024     %           o mode : (0) select a number of Gaussian peaks left unchanged (WARNING : total power less than Pin)
0025     %                    (1) assemble peaks fixing Pmax and N, allowing dN to change
0026     %                    (2) assemble peaks allowing Pmax, N and dN to change
0027     %
0028     %         2) Case where param is empty : to be determined interactively
0029     %
0030     %         3) Case where param is a scalar : param = x + i*y. (default : 2i)
0031     %                 The lobes are selected in the order of decreasing power.
0032     %
0033     %           o case x = 0 : lobes calculated by keeping N and Pmax of
0034     %               each initial selected lobe. dN increased to account for
0035     %               power of neighboring lobes.
0036     %                   - If 0 < y < 1, y is mimimum respective fraction of max power
0037     %                   - If y > 1 (integer), y is number of lobes selected
0038     %                   - If y < -1 (integer), the power cannot be transfered across N=0
0039     %
0040     %           o case x > 0 : lobes calculated by varying N, Pmax and dN to account for
0041     %               neighboring lobes (weighted averages).
0042     %                   - If x < 1, x is mimimum respective fraction of max power
0043     %                   - If x > 1 (integer), x is number of lobes selected
0044     %                   - y (integer >= 0) is the number of completely protected lobes
0045     %
0046     %           o case x < 0 : identical to x = 0, but at the end only |x| (integer)
0047     %               lobes are truly returned. The lobes are picked alternatively in each direction
0048     %             WARNING : in this case the total power is less than Pin
0049     %
0050     %       - opt_optim : use of function lsqnonlin for gaussian lobe fit (default : 0)
0051     %       - sdisplay : display structure, with the following field
0052     %           o inter : interactive mode for selecting parameters
0053     %           o mode : display (0) none (1) text only (2) text and figures (default : 0)
0054     %           o style (optional) : style properties (font...) for interactive windows and figures
0055     %           o opt_gui (optional) : use gui for interactive windows, use as handle if opt_gui is handle
0056     %           o axs (optional) : use axs to plot spectrum (iluke only)
0057     %
0058     %
0059     % OUTPUTS :
0060     %       - yN : values for output lobes
0061     %       - ydN : respective gaussian width
0062     %       - yP : respective power
0063     %       - param_out : processing parameter used to generate spectrum
0064     %       - Pin : initial power used to generate spectrum
0065     %       - zyfac : distribution matrix
0066     %
0067     % by J. Decker and Y. Peysson
0068     %
0069     if nargin < 7,
0070         sdisplay = struct;
0071     end
0072     %
0073     if nargin < 6,
0074         opt_optim = 0;
0075     end
0076     %
0077     if nargin < 5,
0078         param_in = NaN;
0079     end
0080     %
0081     if nargin < 4,
0082         Pin_in = NaN;
0083     end
0084     %
0085     % PARAMETERS
0086     %
0087     if opt_optim == 1 && exist('lsqnonlin') ~= 2,
0088         opt_optim = 0;
0089     end
0090     %
0091     if ~isstruct(param_in) && isnan(param_in),% default values
0092         param_in.dir = 1;
0093         param_in.mode = 0;
0094         param_in.nlobes = 2;
0095     end
0096     %
0097     if isstruct(param_in) && ~isfield(param_in,'dir'),
0098         param_in.dir = 1;
0099     end
0100     %
0101     if isstruct(param_in) && ~isfield(param_in,'mode'),
0102         param_in.mode = 0;
0103     end
0104     %
0105     if isstruct(param_in) && ~isfield(param_in,'nlobes'),
0106         param_in.nlobes = 2;
0107     end
0108     %
0109     % DISPLAY
0110     %
0111     if ~isstruct(sdisplay),
0112         display_mode = sdisplay;
0113         sdisplay = struct;
0114     elseif isfield(sdisplay,'mode'),
0115         display_mode = sdisplay.mode;
0116     else
0117         display_mode = 0;
0118     end
0119     %
0120     if isfield(sdisplay,'inter'),
0121         %
0122         inter_mode = sdisplay.inter;% interactive mode
0123         %
0124     else
0125         %
0126         inter_mode = 1;
0127         %
0128     end     
0129     %
0130     if isfield(sdisplay,'style'),
0131         style = sdisplay.style;
0132     else
0133         style = struct;
0134     end
0135     %
0136     if isfield(sdisplay,'opt_gui'),
0137         %
0138         if ishandle(sdisplay.opt_gui),
0139             selectheight = 175;
0140             panelposition = get(sdisplay.opt_gui,'Position');
0141             %
0142             select.position = [style.panelsep,style.panelsep,panelposition(3) - 2*style.panelsep,selectheight];
0143             info.position = [style.panelsep,2*style.panelsep + selectheight,panelposition(3) - 2*style.panelsep,panelposition(4) - selectheight - 3*style.panelsep];
0144             info1.position = [style.panelsep,2*style.panelsep + selectheight,panelposition(3)/2 - 2*style.panelsep,panelposition(4) - selectheight - 3*style.panelsep];
0145             info2.position = [panelposition(3)/2 + style.panelsep,2*style.panelsep + selectheight,panelposition(3)/2 - 2*style.panelsep,panelposition(4) - selectheight - 3*style.panelsep];
0146             %
0147             select.handle = uipanel('Parent',sdisplay.opt_gui,'Title','','units','pixels','Position',select.position);
0148             info.handle = uicontrol(sdisplay.opt_gui,'Style','text','Position',info.position,'FontName',style.fontname,'FontUnits','pixels','FontSize',style.fontsize,'HorizontalAlignment','left');
0149             info1.handle = uicontrol(sdisplay.opt_gui,'Style','text','Position',info1.position,'FontName',style.fontname,'FontUnits','pixels','FontSize',style.fontsize,'HorizontalAlignment','left');
0150             info2.handle = uicontrol(sdisplay.opt_gui,'Style','text','Position',info2.position,'FontName',style.fontname,'FontUnits','pixels','FontSize',style.fontsize,'HorizontalAlignment','center');
0151             %
0152             opt_gui = select.handle;
0153             opt_info = true;
0154             %
0155             style.textwidth = panelposition(3) - 2*style.panelsep - style.butwidth - 5*style.buthsep;
0156             %
0157         else
0158             opt_gui = sdisplay.opt_gui;
0159             opt_info = false;
0160         end
0161     else
0162         opt_gui = false;
0163         opt_info = false;
0164     end
0165     %
0166     if display_mode >= 1 && inter_mode == 1,
0167         %
0168         infostr = {...
0169         '--------------------------------------------------------',...
0170         '          ALOHA  ->  C3PO parametrization               ',...
0171         '--------------------------------------------------------',...
0172         ' ',...
0173         'Recommended parametrization for aloha -> C3PO :',...
0174         ' ',...      
0175         '  - Multi-lobe model (strong absorption case,... : param.dir = 1, param.mode = 0, param.nlobes >= 2',...
0176         '        Note : adjust param.nlobes to approximately conserve absolute directivity and repartition',...
0177         ' ',...      
0178         '  - Tail-lobes model (weak absorption case,...   : param.dir = 2, param.mode = 1, param.nlobes = 2',...
0179         '        Note : absolute directivity and repartition should be fairly well conserved',...
0180         ' ',...      
0181         '  - Simple model (ideal spectrum, one lobe,...   : param.dir = 1, param.mode = 1, param.nlobes = 1',...           
0182         ' '};     
0183         %
0184         if opt_info,
0185             set(info1.handle,'Visible','off');
0186             set(info2.handle,'Visible','off');
0187             set(info.handle,'String',infostr([5:13]));
0188         else
0189             for iinfo = 1:length(infostr),
0190                 disp(infostr{iinfo});
0191             end
0192         end        
0193         %
0194     end
0195     %
0196     exit_loop = 0;
0197     figs = 0;
0198     %
0199     while exit_loop == 0,
0200         %
0201         zN = zN_in;
0202         zdN = zdN_in;
0203         zdP = zdP_in;
0204         Pin = Pin_in;
0205         param = param_in;
0206         %
0207         if figs == 1,
0208             close(fig1),close(fig2),close(fig3),close(fig4)
0209         end
0210         %
0211         if inter_mode == 1,
0212             param.dir = iselect_jd({'P     + preserve repartition and total power',...
0213                                     'P/N^2 + preserve directivity and repartition',...
0214                                     'P/N^2 + preserve directivity and total power',...
0215                                     'P/N^2 + preserve repartition and total power'},'Sort peaks according to :',opt_gui,style,param.dir + 1) - 1;
0216             %
0217             disp(' ')      
0218         end
0219         %
0220         if isstruct(param),
0221             if param.dir == 0,
0222                 Pin = real(Pin);
0223             elseif param.dir == 1,
0224                 Pin = real(Pin) + 1i;
0225             elseif param.dir == 2,
0226                 Pin = real(Pin) + 2i;
0227             else
0228                 Pin = real(Pin) + 3i;            
0229             end
0230         end
0231         %
0232         zdP = zdP(:).';
0233         zdN = zdN(:).';
0234         zN = zN(:).';
0235         %
0236         if min(zdP) < 0,% correct case of negative power from aloha
0237             disp('warning: min(zdP) < 0, correcting')
0238             disp(' ')      
0239             zdP = zdP - min(zdP);
0240         end
0241         %
0242         global zx zf
0243         %
0244         if length(zdN) == 1 && isnan(zdN), % calculate grid steps if needed
0245             zN_S = [zN(1) - (zN(2)-zN(1))/2, (zN(1:end-1) + zN(2:end))/2, zN(end) + (zN(end)-zN(end-1))/2];
0246             zdN = diff(zN_S);
0247         end
0248         %
0249         % remove edge points with zdP = 0
0250         %
0251         mask = find(zdP > 0,1,'first'):find(zdP > 0,1,'last');
0252         zdP = zdP(mask);
0253         zdN = zdN(mask);
0254         zN = zN(mask);
0255         %
0256         % manages case of successive points with equal power density
0257         %
0258         zptrend = sign(diff(zdP));
0259         %
0260         while any(zptrend == 0),
0261             %
0262             disp('Fixing the case of successive points with equal power density')
0263             %
0264             iflat = find(zptrend == 0);% zppdP(iflat) = zppdP(iflat+1)
0265             %
0266             while ~isempty(iflat),
0267                 if length(iflat) > 1,% could be more than two successive equal values
0268                     imask = find(diff(iflat) > 1,1,'first');% number of successive equal values
0269                     if isempty(imask),% remaining cases are all successive equal values
0270                         imask = length(iflat);
0271                     end
0272                 else
0273                     imask = 1;
0274                 end
0275                 %
0276                 mask = iflat(1):iflat(imask)+1;% indices successive points with equal values
0277                 zN = [zN(1:mask(1)-1),mean(zN(mask)),zN(mask(end)+1:end)];
0278                 zdP = [zdP(1:mask(1)-1),zdP(mask(1)),zdP(mask(end)+1:end)];
0279                 zdN = [zdN(1:mask(1)-1),sum(zdN(mask)),zdN(mask(end)+1:end)];
0280                 %
0281                 iflat(1:imask) = [];
0282                 iflat = iflat - imask;% shift indices
0283                 %
0284             end              
0285             %
0286             zptrend = sign(diff(zdP));
0287             %
0288         end
0289         %
0290         %
0291         %
0292         Ptot = sum(zdP.*zdN);
0293         nz = length(zdP);
0294         %
0295         if ~isnan(real(Pin)),% externally imposed power
0296             zdP = zdP*real(Pin)/Ptot;%set correct total input power (aloha total power calculation may give erratic results)
0297             Ptot = real(Pin);
0298         end
0299         %
0300         if imag(Pin) >= 1,% sort power as P/N^2 instead of P to conserve directivity
0301             zdP = zdP./zN.^2;
0302             Ptot = sum(zdP.*zdN);
0303         end
0304         %
0305         %  display limits
0306         %
0307         Pmindisp = 10^(floor(mean(log10(zdP(find(zdP > 0))))) - 2); 
0308         Pmaxdisp = 10^(ceil(log10(max(zdP))) + 1); 
0309         %
0310         % identify initial lobes
0311         %
0312         zppdP = [0,zdP,0];
0313         %
0314         zptrend = sign(diff(zppdP));
0315         %
0316         if any(zptrend == 0),
0317             error('case of successive points with equal power density not yet treated.')
0318         end
0319         %
0320         zvar = diff(zptrend);
0321         %
0322         ilim = [1,find(zvar == 2),nz]; %indices of lobes separation
0323         %
0324         ny = length(ilim) - 1;% number of initial lobes
0325         %
0326         enforcedir = 0;
0327         %
0328         if isstruct(param),
0329             %
0330             if inter_mode == 1,
0331                 param.mode = iselect_jd({'select a number of Gaussian peaks left unchanged (WARNING : total power less than Pin)',...
0332                                          'assemble peaks fixing Pmax and N, allowing dN to change',...
0333                                          'assemble peaks allowing Pmax, N and dN to change'},'Mode for peak processing :',opt_gui,style,param.mode + 1) - 1;
0334                 disp(' ')      
0335                 %
0336                 param.nlobes = iselect_jd(num2cell(num2str((1:ny).'),2).','Number of lobes to be selected :',opt_gui,style,param.nlobes,'popupmenu');
0337                 disp(' ')      
0338             end
0339             %
0340             if param.mode == 0,
0341                 param = - param.nlobes - Inf*1i;
0342             elseif param.mode == 1,
0343                 %
0344                 if param.nlobes == 1,
0345                     param = - 1 - 2i;
0346                 else
0347                     param = - param.nlobes*1i;
0348                 end
0349                 %
0350             else
0351                 if param.nlobes > 1,
0352                     enforcedir = 1;
0353                 end
0354                 param = param.nlobes;
0355             end
0356         end
0357         %
0358         if imag(param) <= -2,% do not allow power to be transfered through N=0, usual mode
0359             enforcedir = 1;
0360             param = real(param) - 1i*imag(param);
0361         end
0362         %
0363         if real(param) <= 0,% Peak position and height are maintained, dNpar is adjusted
0364             if imag(param) < 1,
0365                 Pfacmin = imag(param);
0366                 nlobesmax = ny;
0367             else
0368                 Pfacmin = 0;
0369                 nlobesmax = imag(param);
0370             end        
0371             %
0372             nlobespro = 0;
0373             opt_peak = 0;
0374             %
0375         else % Peak position and height are varied for non-protected peaks
0376             if real(param) < 1,
0377                 Pfacmin = real(param);
0378                 nlobesmax = ny;
0379             else
0380                 Pfacmin = 0;
0381                 nlobesmax = real(param);
0382             end
0383             %
0384             % a number of (main) peaks can be protected (they keep the original N, P and dN)
0385             %
0386             % they are not accounted for in nlobesmax (nlobes = nlobesmax + nlobespro).
0387             % if Pfacmin is used, the fraction is calculated on the remaining power
0388             %
0389             nlobespro = fix(imag(param));
0390             opt_peak = 1;
0391             %
0392         end
0393         %
0394         zydP = repmat(zdP.',[1,ny]);
0395         zydN = repmat(zdN.',[1,ny]);
0396         zyN = repmat(zN.',[1,ny]);
0397         %
0398         zyind = (1:nz).'*ones(1,ny);
0399         zyilimm = repmat(ilim(1:end-1),[nz,1]);
0400         zyilimp = repmat(ilim(2:end),[nz,1]);
0401         %
0402         zyfac = zeros(nz,ny);% distribution matrix
0403         zyfac(zyind > zyilimm & zyind < zyilimp) = 1;
0404         zyfac(zyind == zyilimm | zyind == zyilimp) = 1/2;%half power belongs to lobe at separation points
0405         zyfac(zyind == 1 | zyind == nz) = 2*zyfac(zyind == 1 | zyind == nz);%full power belongs to lobe at end points
0406         %
0407         yP = sum(zydP.*zydN.*zyfac);
0408         %
0409         % initial lobes
0410         %
0411         if opt_peak == 1,% recalculate N//
0412             %
0413             yN = sum(zydP.*zydN.*zyfac.*zyN)./yP;
0414             ydN = sqrt(2*sum(zydP.*zydN.*zyfac.*(zyN - repmat(yN,[nz,1])).^2)./yP);
0415             if opt_optim,
0416                 for iy = 1:ny,
0417                     zx = zN - yN(iy);
0418                     zf = zdP.*zyfac(:,iy).'/yP(iy);
0419                     ydN(iy) = lsqnonlin(@gauss_squares_jd,ydN(iy));
0420                 end
0421             end
0422             %
0423         else% keep N//,
0424             %
0425             ipeak = find(zvar == -2);
0426             yN = zN(ipeak);
0427             ydN = yP/sqrt(pi)./zdP(ipeak);% adjust dN to account for peak assembly
0428             %
0429         end
0430         %
0431         if abs(sum(yP) - Ptot)/Ptot > eps('double')*nz,
0432             error('Power was not conserved')
0433         end
0434         %
0435         if nlobespro > 0 && ny > nlobespro,% keep a number of protected peaks
0436             %
0437             [yPsort,iysort] = sort(yP);
0438             %
0439             iypro = iysort(ny + 1 - nlobespro:ny);
0440             %
0441             yNpro = yN(iypro);
0442             ydNpro = ydN(iypro);
0443             yPpro = yP(iypro);
0444             zyfacpro = zyfac(:,iypro);
0445             %
0446             yN(iypro) = [];
0447             ydN(iypro) = [];
0448             yP(iypro) = [];
0449             zyfac(:,iypro) = [];
0450             %
0451             Ptot = Ptot - sum(yPpro);
0452             %
0453             ny = ny - nlobespro;
0454             %
0455         else
0456             %
0457             yNpro = yN([]);
0458             ydNpro = ydN([]);
0459             yPpro = yP([]);
0460             zyfacpro = zyfac(:,[]);
0461             %
0462         end
0463         %
0464         yP0 = yP;
0465         %
0466         % loop to agregate the power of minor peaks
0467         %
0468         while ny > nlobesmax || min(yP)/Ptot < Pfacmin,
0469             %
0470             iy = find(yP == min(yP),1,'first');
0471             %
0472             if iy == 1,
0473                 %
0474                 iyp = 2;
0475                 %
0476                 if enforcedir == 1 && yN(iyp)*yN(iy) < 0,% forbid changes in directivity
0477                     %
0478                     yNpro = [yNpro,yN(iy)];
0479                     ydNpro = [ydNpro,ydN(iy)];
0480                     yPpro = [yPpro,yP(iy)];
0481                     zyfacpro = [zyfacpro,zyfac(:,iy)];
0482                     %
0483                     if opt_peak == 0,% recalculate dN// to keep Pmax.
0484                         %
0485                         ydNpro(end) = ydN(iy).*yP(iy)./yP0(iy);
0486                         %
0487                     end
0488                     %
0489                     Ptot = Ptot - yP(iy);
0490                     %
0491                     yN(iy) = [];
0492                     ydN(iy) = [];
0493                     yP(iy) = [];
0494                     yP0(iy) = [];
0495                     zyfac(:,iy) = [];
0496                     %
0497                     ny = ny - 1;
0498                     %
0499                     nlobesmax = nlobesmax - 1;
0500                     nlobespro = nlobespro + 1;
0501                     %
0502                     continue
0503                     %
0504                 end
0505                 %
0506                 zyfac(:,iyp) = zyfac(:,iyp) + zyfac(:,iy);
0507                 yP(iyp) = sum(zdP.*zdN.*zyfac(:,iyp).');
0508                 %
0509                 if opt_peak == 1,% recalculate N//
0510                     yN(iyp) = sum(zdP.*zdN.*zyfac(:,iyp).'.*zN)./yP(iyp);
0511                     ydN(iyp) = sqrt(2*sum(zdP.*zdN.*zyfac(:,iyp).'.*(zN - yN(iyp)).^2)/yP(iyp));
0512                     if opt_optim,
0513                         zx = zN - yN(iyp);
0514                         zf = zdP.*zyfac(:,iyp).'/yP(iyp);
0515                         ydN(iyp) = lsqnonlin(@gauss_squares_jd,ydN(iyp));
0516                     end
0517                 end
0518                 %
0519             elseif iy == ny,
0520                 %
0521                 iym = ny - 1;
0522                 %
0523                 if enforcedir == 1 && yN(iym)*yN(iy) < 0,% forbid changes in directivity
0524                     %
0525                     yNpro = [yNpro,yN(iy)];
0526                     ydNpro = [ydNpro,ydN(iy)];
0527                     yPpro = [yPpro,yP(iy)];
0528                     zyfacpro = [zyfacpro,zyfac(:,iy)];
0529                     %
0530                     if opt_peak == 0,% recalculate dN// to keep Pmax.
0531                         %
0532                         ydNpro(end) = ydN(iy).*yP(iy)./yP0(iy);
0533                         %
0534                     end
0535                     %
0536                     Ptot = Ptot - yP(iy);
0537                     %
0538                     yN(iy) = [];
0539                     ydN(iy) = [];
0540                     yP(iy) = [];
0541                     yP0(iy) = [];
0542                     zyfac(:,iy) = [];
0543                     %
0544                     ny = ny - 1;
0545                     %
0546                     nlobesmax = nlobesmax - 1;
0547                     nlobespro = nlobespro + 1;
0548                     %
0549                     continue
0550                     %
0551                 end
0552                 %
0553                 zyfac(:,iym) = zyfac(:,iym) + zyfac(:,iy);
0554                 yP(iym) = sum(zdP.*zdN.*zyfac(:,iym).');
0555                 %
0556                 if opt_peak == 1,% recalculate N//
0557                     yN(iym) = sum(zdP.*zdN.*zyfac(:,iym).'.*zN)./yP(iym);
0558                     ydN(iym) = sqrt(2*sum(zdP.*zdN.*zyfac(:,iym).'.*(zN - yN(iym)).^2)/yP(iym));
0559                     if opt_optim,
0560                         zx = zN - yN(iym);
0561                         zf = zdP.*zyfac(:,iym).'/yP(iym);
0562                         ydN(iym) = lsqnonlin(@gauss_squares_jd,ydN(iym));
0563                     end
0564                 end
0565                 %
0566             else
0567                 %
0568                 iym = iy - 1;
0569                 iyp = iy + 1;
0570                 %
0571                 if enforcedir == 1 && yN(iym)*yN(iyp) < 0,% forbid changes in directivity
0572                     %
0573                     Nlim = 0;
0574                     %
0575                 else
0576                     %
0577                     Nlim = (yN(iym) + yN(iyp))/2;
0578                     %
0579                 end
0580                 %
0581                 izlimm = find(zN_S < Nlim,1,'last');
0582                 %
0583                 zyfac(1:izlimm-1,iym) = zyfac(1:izlimm-1,iym) + zyfac(1:izlimm-1,iy);
0584                 zyfac(izlimm+1:end,iyp) = zyfac(izlimm+1:end,iyp) + zyfac(izlimm+1:end,iy);
0585                 %
0586                 zyfac(izlimm,iym) = zyfac(izlimm,iym) + zyfac(izlimm,iy)*(Nlim - zN_S(izlimm))/zdN(izlimm);
0587                 zyfac(izlimm,iyp) = zyfac(izlimm,iyp) + zyfac(izlimm,iy)*(zN_S(izlimm + 1) - Nlim)/zdN(izlimm);
0588                 %
0589                 yP(iym) = sum(zdP.*zdN.*zyfac(:,iym).');
0590                 yP(iyp) = sum(zdP.*zdN.*zyfac(:,iyp).');
0591                 %
0592                 if opt_peak == 1,% recalculate N//
0593                     yN(iym) = sum(zdP.*zdN.*zyfac(:,iym).'.*zN)./yP(iym);
0594                     ydN(iym) = sqrt(2*sum(zdP.*zdN.*zyfac(:,iym).'.*(zN - yN(iym)).^2)/yP(iym));
0595                     if opt_optim,
0596                         zx = zN - yN(iym);
0597                         zf = zdP.*zyfac(:,iym).'/yP(iym);
0598                         ydN(iym) = lsqnonlin(@gauss_squares_jd,ydN(iym));
0599                     end
0600                     %
0601                     yN(iyp) = sum(zdP.*zdN.*zyfac(:,iyp).'.*zN)./yP(iyp);
0602                     ydN(iyp) = sqrt(2*sum(zdP.*zdN.*zyfac(:,iyp).'.*(zN - yN(iyp)).^2)/yP(iyp));
0603                     if opt_optim,
0604                         zx = zN - yN(iyp);
0605                         zf = zdP.*zyfac(:,iyp).'/yP(iyp);
0606                         ydN(iyp) = lsqnonlin(@gauss_squares_jd,ydN(iyp));
0607                     end
0608                 end
0609                 %
0610             end
0611             %
0612             yP0(iy) = [];
0613             yP(iy) = [];
0614             yN(iy) = [];
0615             ydN(iy) = [];
0616             zyfac(:,iy) = [];
0617             %
0618             ny = ny - 1;    
0619             %
0620         end
0621         %
0622         if opt_peak == 0,% recalculate dN// to keep Pmax.
0623             %
0624             ydN = ydN.*yP./yP0;
0625             %
0626         end
0627         %
0628         % add protected peaks
0629         %
0630         yN = [yNpro,yN];
0631         ydN = [ydNpro,ydN];
0632         yP = [yPpro,yP];
0633         zyfac = [zyfacpro,zyfac];
0634         %
0635         ny = ny + nlobespro;
0636         %
0637         % select number of peaks if real(param) < 0
0638         %
0639         if real(param) < 0,
0640             %
0641             ny = min([ny,-real(param)]);
0642             %
0643             [dummy,iy] = sort(yP.^2,2,'descend');% keep this to choose with Pin whether to work with directivity
0644             %
0645             iys = iy(1);
0646             iy(1) = [];
0647             %
0648             for iiy = 2:ny,
0649                 j = 1;
0650                 while j <= length(iy),
0651                     if yN(iy(j))*yN(iys(end)) < 0,% pick alternatively in each direction
0652                         iys(iiy) = iy(j);
0653                         iy(j) = [];
0654                         break
0655                     else
0656                         j = j + 1;
0657                     end
0658                 end
0659                 %
0660                 if j > length(iy),% one direction is exhausted
0661                     iys(iiy) = iy(1);
0662                     iy(1) = [];
0663                 end
0664             end
0665             %
0666             yP = yP(iys);
0667             yN = yN(iys);
0668             ydN = ydN(iys);
0669             zyfac = zyfac(:,iys);
0670             %
0671         end
0672         %
0673         % sort output by N//
0674         %
0675         [yN,iy] = sort(yN);
0676         ydN = ydN(iy);
0677         yP = yP(iy);
0678         zyfac = zyfac(:,iy);
0679         %
0680         if imag(Pin) >= 1,% readjust power
0681             yP = yP.*yN.^2;
0682             zdP = zdP.*zN.^2;
0683             Ptot = sum(zdP.*zdN);
0684         end    
0685         %
0686         if imag(Pin) == 2,% readjust power keeping directivity and total power
0687             %
0688             maskp = yN > 0;
0689             maskm = yN < 0;
0690             %
0691             Ptotp = sum(yP(maskp));
0692             Ptotm = sum(yP(maskm));
0693             %
0694             Pdirp = sum(yP(maskp)./yN(maskp).^2);
0695             Pdirm = sum(yP(maskm)./yN(maskm).^2);
0696             %
0697             facp = (Ptot*Pdirm + Ptotm*(Pdirp - Pdirm))/(Pdirm*Ptotp + Pdirp*Ptotm);
0698             facm = (Ptot*Pdirp + Ptotp*(Pdirm - Pdirp))/(Pdirm*Ptotp + Pdirp*Ptotm);
0699             %
0700             yP(maskp) = yP(maskp)*facp;
0701             yP(maskm) = yP(maskm)*facm; 
0702             %
0703     %         Ptotp*facp + Ptotm*facm = Ptot
0704     %         Pdirp*facp - Pdirm*facm = (Pdirp*Ptot*Pdirm + Pdirp*Ptotm*(Pdirp - Pdirm) - Pdirm*Ptot*Pdirp - Pdirm*Ptotp*(Pdirm - Pdirp))/(Pdirm*Ptotp + Pdirp*Ptotm)
0705     %                                 = Pdirp - Pdirm
0706             %
0707             if opt_peak == 0 && real(param) == 0,% recalculate dN// to keep Pmax.
0708                 %
0709                 ydN(maskp) = ydN(maskp)*facp;
0710                 ydN(maskm) = ydN(maskm)*facm;
0711                 %
0712             end
0713         end
0714         %
0715         if imag(Pin) == 3,% readjust power keeping total power
0716             %
0717             fac = Ptot/sum(yP);
0718             yP = yP*fac;
0719             %
0720             if opt_peak == 0 && real(param) == 0,% recalculate dN// to keep Pmax.
0721                 %
0722                 ydN = ydN*fac;
0723                 %
0724             end
0725             %
0726         end
0727         %
0728         % display results
0729         %
0730         if display_mode >= 1,
0731             %
0732             Npeak = zN(find(zdP == max(zdP),1,'first'));
0733             %
0734             zmaskp = zN > 0;
0735             zmaskm = zN < 0;
0736             %
0737             ymaskp = yN > 0;
0738             ymaskm = yN < 0;
0739             %
0740             % Ideal weighted directivity (RC = 0), X. Litaudon and D. Moreau 1990 Nucl. Fusion 30 471
0741             %
0742             di = Npeak^2*(sum(zdP(zmaskp)./zN(zmaskp).^2.*zdN(zmaskp)) - sum(zdP(zmaskm)./zN(zmaskm).^2.*zdN(zmaskm)))/Ptot;
0743             df = Npeak^2*(sum(yP(ymaskp)./yN(ymaskp).^2) - sum(yP(ymaskm)./yN(ymaskm).^2))/Ptot;
0744             %
0745             pti = Ptot;
0746             ptf = sum(yP);
0747             %
0748             dri = di;
0749             drf = df*pti/ptf;
0750             %
0751             %
0752             % Module directivity (RC = 0), J. Hillairet et al. Nucl. Fusion 50 (2010) 125010
0753             %
0754             ri = sum(zdP(zmaskp).*zdN(zmaskp))/pti;
0755             rf = sum(yP(ymaskp))/ptf; 
0756             %
0757             infostr = {...
0758             ['Ideal weighted current drive directivity (RC = 0), X. Litaudon and D. Moreau 1990 Nucl. Fusion 30 471'],...
0759             ['------------------------------------------------------------------------------------------------'],...
0760             ['Initial   directivity : ',num2str(di)],...
0761             ['Processed directivity : ',num2str(df)],...        
0762             ' ',...
0763             ['Initial   power (W) : ',num2str(pti)],...
0764             ['Processed power (W) : ',num2str(ptf)],...        
0765             ' ',...
0766             ['Relative directivity'],...
0767             ['------------------------------------------------------------------------------------------------'],...
0768             ['Initial   rel. directivity : ',num2str(dri)],...
0769             ['Processed rel. directivity : ',num2str(drf)],...        
0770             ' ',...
0771             ['Power directivity, J. Hillairet et al. Nucl. Fusion 50 (2010) 125010'],...
0772             ['------------------------------------------------------------------------------------------------'],...
0773             ['Initial   repartition : ',num2str(ri)],...
0774             ['Processed repartition : ',num2str(rf)],...        
0775             ' ',...
0776             '   yN     ydN     yP   ',...
0777             '-----------------------',...
0778             num2str([yN.',ydN.',yP.'],3),...
0779             ' '};
0780             %
0781             if opt_info,
0782                 set(info.handle,'Visible','off');
0783                 set(info1.handle,'Visible','on');
0784                 set(info2.handle,'Visible','on');
0785                 set(info1.handle,'String',infostr([1,2,3,4, 6,7, 9,10,11,12,14,15,16,17]));
0786                 set(info2.handle,'String',infostr(13:end-1));
0787             else
0788                 for iinfo = 1:length(infostr),
0789                     disp(infostr{iinfo});
0790                 end
0791             end        
0792             %
0793         end
0794         %
0795         if display_mode >= 2,
0796             %
0797             dx = min(ydN)/10;
0798             x = (min(yN) - 4*max(ydN)):dx:(max(yN) + 4*max(ydN));
0799             nx = length(x);
0800             %
0801             xyN = repmat(yN,[nx,1]);
0802             xyP = repmat(yP,[nx,1]);
0803             xydN = repmat(ydN,[nx,1]);
0804             xy = repmat(x.',[1,ny]);
0805             %
0806             xyf = xyP.*exp(-(xy - xyN).^2./xydN.^2)./sqrt(pi)./xydN;
0807             xf = sum(xyf,2);
0808             %
0809             ylim = [Pmindisp,Pmaxdisp];
0810             %
0811             if real(Pin) == 1,
0812                 fac = 1;
0813                 ylab = 'dP/dn_z (norm.)';
0814             else
0815                 fac = 1e6;
0816                 ylab = 'dP/dn_z (MW)';
0817             end            
0818             %
0819             if isfield(sdisplay,'axs'),
0820                 ax = sdisplay.axs(1);
0821                 set(ax,'Visible','on');
0822                 red = 1;
0823                 lspace = NaN;
0824                 bspace = NaN;
0825             else
0826                 fig1 = figure('Name','Gaussian fit vs. ALOHA');
0827                 ax = gca;
0828                 red = 0.9;
0829                 lspace = 0.7;
0830                 bspace = 0.5;
0831             end
0832             %
0833             leg = {'ALOHA','Gaussian fit','Location','NorthWest'};
0834             graph1D_jd(zN,zdP/fac,0,0,'n_z',ylab,'',NaN,'',NaN,'-','none','b',2,style,ax,red,lspace,bspace);
0835             graph1D_jd(x,xf/fac,0,0,'','','',leg,'',NaN,'--','none','r',2,style,ax);
0836             %
0837             if ~isfield(sdisplay,'axs'),
0838                 fig2 = figure('Name','Gaussian fit vs. ALOHA');
0839                 graph1D_jd(zN,zdP/fac,0,1,'n_z',ylab,'',NaN,'',ylim/fac,'-','none','b',2,style,gca,red,lspace,bspace);
0840                 graph1D_jd(x,xf/fac,0,1,'','','',leg,'',ylim/fac,'--','none','r',2);
0841                 %
0842                 fig3 = figure('Name','Gaussian fit calculation');
0843                 leg = {'Gaussian','Fit'};
0844                 graph1D_jd(xy(:),xyf(:)/fac,0,0,'n_z',ylab,'',NaN,'',NaN,'-','none','b',2,style,gca,red,lspace,bspace);
0845                 graph1D_jd(x,xf/fac,0,0,'','','',leg,'',NaN,'--','none','r',2);
0846                 %
0847                 fig4 = figure('Name','Gaussian fit calculation');
0848                 graph1D_jd(xy(:),xyf(:)/fac,0,1,'n_z',ylab,'',NaN,'',ylim/fac,'-','none','b',2,style,gca,red,lspace,bspace);
0849                 graph1D_jd(x,xf/fac,0,1,'','','',leg,'',ylim/fac,'--','none','r',2);
0850                 %
0851                 figs = 1;
0852                 %
0853             end
0854             %
0855         end
0856         %
0857         if inter_mode == 0,
0858             exit_loop = 1;
0859         elseif isfield(sdisplay,'exit_loop') && sdisplay.exit_loop,
0860             exit_loop = sdisplay.exit_loop;
0861         else
0862             exit_loop = 2 - iselect_jd({'Accept this process','Start again'},'Do you want to :',opt_gui,style,1);
0863             disp(' ')      
0864         end
0865         %
0866         if exit_loop == 0 && isfield(sdisplay,'axs'),% start again - clear new plot
0867             cla;% clear graph of lines
0868             chs = get(get(ax,'parent'),'children');% list of panel elements
0869             delete(chs(1));% delete legend
0870         end
0871         %
0872         if exit_loop == 0 && opt_info,% start again - clear disp
0873             set(info.handle,'Visible','on');
0874             set(info1.handle,'Visible','off');
0875             set(info2.handle,'Visible','off');
0876         end
0877     end
0878     %
0879     if opt_info,
0880         delete(select.handle);
0881         delete(info.handle);
0882         delete(info1.handle);
0883         delete(info2.handle);
0884     end
0885     %
0886     if imag(param) >= 2 && enforcedir == 1;
0887         param = real(param) - 1i*imag(param);
0888     end 
0889     %
0890     param_out.calc = param;
0891     param_out.dir = imag(Pin);
0892     %
0893     if imag(param) == -Inf && real(param) <= -1,
0894         param_out.mode = 0;
0895         param_out.nlobes = -real(param);
0896     elseif param == - 1 - 2i;
0897         param_out.mode = 1;
0898         param_out.nlobes = 1;
0899     elseif imag(param) < -1 && real(param) == 0,
0900         param_out.mode = 1;
0901         param_out.nlobes = -imag(param);
0902     elseif imag(param) == 0 && real(param) >= 1,
0903         param_out.mode = 2;
0904         param_out.nlobes = real(param);
0905     else
0906         param_out.mode = NaN;
0907         param_out.nlobes = NaN;
0908     end
0909     %
0910 end       
0911 %
0912 function [zdf] = gauss_squares_jd(dN)
0913     %
0914     global zx zf
0915     %
0916     zdf = zf - exp(-zx.^2/dN^2)/sqrt(pi)/dN;
0917     %
0918 end
0919 
0920 
0921 
0922 
0923 
0924 
0925