Generating ASCOT5 input entirely from TRANSP about ascot5 HOT 7 OPEN

  function a=transp2ascot4distribution(transpCDF,ascotH5)
%
% We cheat by putting all the variuous RZ-points into a single column in 
%
%
   distName='rzPitchEdist';

   disp('Reading transp...')
   tra=read_transp_netcdf_aug(transpCDF);
   
   if(mod(numel(tra.grid.r),2) ~= 0)
      error('Our scheme to have two R columns and many lines failed!') 
   end
   
   disp('Converting...')
   d.version    = 3;
   d.ndims      = 6;
   d.nslots     = round([2 numel(tra.grid.r)/2 numel(tra.pitch) numel(tra.energy) 1 1]);
   d.vectlength = 1;
   
   ordinate = permute(tra.dist,[3 2 1]);
   d.ordinate(1,1,:,:,:) = ordinate(1:(numel(tra.grid.r)/2),    :,:);
   d.ordinate(1,2,:,:,:) = ordinate((numel(tra.grid.r)/2+1):end,:,:);
    
   d.abscissa_edges={               ...
       [1.0 2.0 3.0],                   ...
       linspace(1,2,d.nslots(2)+1), ...
       tra.pitch_boundary,          ...
       tra.energy_boundary,         ...
       [0.0 1.0],                   ...
       [0.5 1.5]                    };
   d.abscissa_units =  {'m',  'm',  'vpar/vtot',  'eV',     's',    ''};
   d.abscissa_names =  {'R',  'z',  'pitch',      'energy', 'time', 'species'};
   d.ordinate_name  =  {'density'};
   d.ordinate_unit  =  {'m^{-3}eV^{-1}'};

   
   a.species.testParticle.anum   = 2;
   a.species.testParticle.znum   = 1;
   a.species.testParticle.mass   = 1.6605e-27;
   a.species.testParticle.charge = 1.6022e-19;

   
   

   if(nargout > 0)
       a.dists.(distName) = d;
   end
   
   if(nargin > 1)
        disp('Writing...')
        write_ascot4_distribution(d,distName,ascotH5);
        write_ascot4_species(a.species.testParticle, ascotH5);
   end

end

function tra=read_transp_netcdf_aug(filename)
% This file reads a nubeam file, such as: 29226A02_fi_1.cdf


tra.filemame=filename;
tra.grid.z    =ncread(filename,'Z2D')/100;
tra.grid.r    =ncread(filename,'R2D')/100;
tra.grid.theta=ncread(filename,'TH2D');
tra.grid.X    =ncread(filename,'X2D');
tra.grid.volumes =ncread(filename,'BMVOL');

tra.pitch     =ncread(filename,'A_D_NBI');
tra.pitch_boundary=ncread(filename,'AB_D_NBI'); 
tra.energy    =ncread(filename,'E_D_NBI'); %eV
tra.energy_boundary=ncread(filename,'EB_D_NBI'); 
tra.energy_units='eV';

tra.dist      =ncread(filename,'F_D_NBI') * 0.5 ;
% Caveat about the normalisation of FBM w.r.t. the pitch angle
%The TRANSP coordinate for the distribution function FBM is dA = dOmega/4pi as reported in the netCDF file label.
%Since there is no FBM dependence on theta, the integral over all pitch angles is
%2pi/4pi \int_0^2pi FBM*sin(theta) d(theta) 
%which is equivalent to consider FBM as the probability density with respect to the coordinate g = cos(theta)/2, i.e.
%df/dg = FBM
%Other codes use p = cos theta = 2*g as direct pitch angle coordinate. In that case their pitch angle probability density is:
%df/dp = df/d(2g) = 0.5 * df/dg = 0.5 * FBM 
% https://www.aug.ipp.mpg.de/aug/manuals/transp/fbm/pitch.html

% d(cos(theta))=sin(theta)d(theta)


tra.dist_units=ncreadatt(filename,'F_D_NBI','units');

tra.time=ncread(filename,'TIME');

tra.runid=ncread(filename,'TRANSP_RUNID');
disp(['Transp runid: "' tra.runid' '"'])


tra.n_fi=squeeze(sum(sum(tra.dist,2),1));
tra.n_fi=tra.n_fi*100^3; % 1/cm^3 => 1/m^3
tra.n_fi=tra.n_fi*(tra.pitch_boundary(2)-tra.pitch_boundary(1));
tra.n_fi=tra.n_fi*(tra.energy_boundary(2)-tra.energy_boundary(1));



r=linspace(min(tra.grid.r),max(tra.grid.r),30);
z=linspace(min(tra.grid.z),max(tra.grid.z),50);

d = griddata(tra.grid.r,tra.grid.z,tra.n_fi,r,z');
set(pcolor(r,z,d),'linestyle','none');
hold on;
scatter(tra.grid.r,tra.grid.z,50,tra.n_fi,'filled');
axis image
box on
xlabel('R (m)')
ylabel('z (m)')
cbh=colorbar;
%ylabel(cbh,'1/(m^3 d(omega/4pi))')
ylabel(cbh,'1/(m^3)')
title('Fast ion density')

colormap(jetW)

tra.volumeAverage = sum(tra.n_fi .* tra.grid.volumes) ./ sum(tra.grid.volumes);
fprintf('Volume averaged density is %e 1/m^3\n',tra.volumeAverage)

tra.vdist = sum(tra.dist.*permute(repmat(tra.grid.volumes,[1,75,50]),[2 3 1]),3);
end

function [a,tra]=transp_4dPE_dist_2_ascot_4dPE_dist(transpCDF,ascotH5)

    %% PARAMETERS  
    
    constants
  
    
    nR       = 30;
    nZ       = 42;
  
    
    mass   = m_d;
    charge = e;
    anum   = 2;
    znum   = 1;
    origin = 1;
    
    %% DISTRIBUTION

   distName='rzPitchEdist';

   disp('Reading transp...')
   tra=read_transp_netcdf_aug(transpCDF);
   drawnow   
   
   disp('Initializing...')
   
   
    nE     = numel(tra.energy);
    nP     = numel(tra.pitch);
    nT       = 1;
    nSpecies = 1;
    nOrdDim  = 1;

     
   tra.grid.rmin=min(tra.grid.r);
   tra.grid.rmax=max(tra.grid.r);
   rmin =  tra.grid.rmin - 0.1*(tra.grid.rmax - tra.grid.rmin);
   rmax =  tra.grid.rmax + 0.1*(tra.grid.rmax - tra.grid.rmin);

   tra.grid.zmin=min(tra.grid.z);
   tra.grid.zmax=max(tra.grid.z);
   zmin =  tra.grid.zmin - 0.1*(tra.grid.zmax - tra.grid.zmin);
   zmax =  tra.grid.zmax + 0.1*(tra.grid.zmax - tra.grid.zmin);
   
   
   d.version    = 3;
   d.ndims      = 6;
   d.nslots     = round([nR nZ nP nE  nT nSpecies]');
   d.vectlength = nOrdDim;
   
    
   d.abscissa_edges={               ...
       linspace( rmin,rmax,d.nslots(1)+1),               ...
       linspace( zmin,zmax,d.nslots(2)+1),               ...
       tra.pitch_boundary,                               ...
       tra.energy_boundary * e,                          ...
       [0.0 1.0],                   ...
       [0.5 1.5]                    };
   d.abscissa_units =  {'m',  'm',  '',  'J',     's',        ''};
   d.abscissa_names =  {'R',  'z',  'pitch',  'energy',  'time', 'species'};
   d.ordinate_name  =  {'density'};
   d.ordinate_unit  =  {'s^3/m^6'};

   for iabs=1:d.ndims
    d.abscissa{iabs}=mean(...
        [ d.abscissa_edges{iabs}(1:(end-1)) ... 
          d.abscissa_edges{iabs}(2:end) ],2);
    d.abscissa_size{iabs}=...   
        diff(d.abscissa_edges{iabs});            
   end
   
   
   disp('Precalculating r,z weights for each transp slot for each ascot rz-gridpoint')
   
   r  = 0.5* ( d.abscissa_edges{1}(1:(end-1)) + d.abscissa_edges{1}(2:(end)));
   z  = 0.5* ( d.abscissa_edges{2}(1:(end-1)) + d.abscissa_edges{2}(2:(end)));
   
   nRZ=numel(tra.grid.r);
   w=zeros(nRZ,nR,nZ);
   for iRZ=1:nRZ;
       values=zeros(size(tra.grid.r));
       values(iRZ)=1.0;
       w(iRZ,:,:) = griddata(tra.grid.r,tra.grid.z,values,r,z')';
       if(false)
           clf
           pcolor(r,z,squeeze(w(iRZ,:,:))')
           hold on
           plot(tra.grid.r(iRZ),tra.grid.z(iRZ),'rx')
           caxis([0 1])
           colormap(jetW)
           drawnow
           pause(1)
       end
   end
   
   w(isnan(w)) = 0.0;
   
   for iR=nR:-1:1
       for iZ=nZ:-1:1
           rz_ind{iR,iZ}=find(w(:,iR,iZ) ~= 0);
           rz_wei{iR,iZ}=w( rz_ind{iR,iZ} ,iR,iZ);
       end
   end
   
   if(false)
       for iR=1:2:nR
           for iZ=1:7:nZ
               clf
               scatter(tra.grid.r,tra.grid.z,20,w(:,iR,iZ),'filled');
               hold on
               plot(r(iR),z(iZ),'rx');
               caxis([0 1])
               %           colormap(jetW)
               title(sum(w(:,iR,iZ)))
               drawnow          
               pause(1)
           end
       end
   end

   disp('Converting transp to ASCOT4 Pitch-Energy...')

   
   
   
   
   d.ordinate = zeros(d.nslots');
   
   
   
    for iZ=1:nZ
       fprintf('.');
        for iR=1:nR
%            for iE=1:nE
%                for iP=1:nP

            if(~isempty(rz_wei{iR,iZ}))
                rzw=repmat(reshape(rz_wei{iR,iZ},1,1,[]),[nE nP 1]);
                d.ordinate(    iR,iZ,:,:,1,1) = ...
                        sum( (tra.dist(:,:,rz_ind{iR,iZ})  ...
                        .* rzw),3 )' ;
            end
%                end
%            end
        end
    end
   fprintf('\n');

   
   
     
   
   
   %% Units and Jacobian
   
   J=1/(100^-3 .* e) ; % 1/(cm^3 eV) => 1/(m^3 eV)
   %keyboard
   

   
   
   d.ordinate = reshape(d.ordinate.* J,[1 size(d.ordinate) ]);
   
  
   
   
  
   
   
   %% TEST PARTICLE SPECIES   
   
   a.species.testParticle.anum   = anum;
   a.species.testParticle.znum   = znum;
   a.species.testParticle.mass   = mass;
   a.species.testParticle.charge = charge;
   a.species.testParticle.origin = origin;

   %% Also create vpa,vpe distribution
   
   
   a.dists.(distName) = d;
   disp('Converting ASCOT4 Pitch-Energy to ASCOT4 Vpa-Vpe...')

   so=convert_ascot4_4dDist_PE_2_VpaVpe(a);
   %so.dists
   
   %% OUTPUT VARIABLES

   if(nargout > 0)
       a.dists.rzVDist = so.dists.rzVDist;
   end
   
   if(nargin > 1)
        disp(['Writing... to file "' ascotH5  '".'])
        write_ascot4_distribution(d,distName,ascotH5);
        write_ascot4_distribution(so.dists.rzVDist,'rzVDist',ascotH5);
        write_ascot4_species(a.species.testParticle, ascotH5);
   end

   %   keyboard


   dv2=diff(d.abscissa_edges{3}(1:2)) *... 
       diff(d.abscissa_edges{4}(1:2));
   
   figure()
   plot2dgrid(d.abscissa_edges{1},d.abscissa_edges{2}, ...
              dv2*squeeze(sum(sum(d.ordinate,5),4)));
   colormap(jetW)
   colorbar
   axis image


end


function [w,ind]=peWeightsInd(x0,xvec)
  
    if(x0 < xvec(1))
        ind= 1;
        w  = 1;
        return
    end
    if(x0 > xvec(end))
        ind= numel(xvec);
        w  = 1;
        return
    end
        
    
    [~,in]=min(abs(xvec-x0));
    
    xn=xvec(in);
    
    if( xn == x0)
      ind = in;
      w   = 1.0;
      return
    end

    if(x0 > xn)
      ind = [ in in+1 ] ;
    else
      ind = [ in-1 in ] ;        
    end
    
    if(min(ind) == 0)
      keyboard  
    end
    
    w(2)  = (x0 - xvec(ind(1))) / ( xvec(ind(2)) - xvec(ind(1)))  ;
    w(1)  = 1.0 - w(2);
    
    
end