function [v,vx,vy,vz,vxx,vxy,vxz,vyy,vyz,vzz,slay,sss,lonv,latv] = newton_trig_int_ver1(lon,lat,density,R1,R2,dSp,G,fntl,ngl,lays,nres,Rh,tlon,tlat,tg,ief,iefs,layerss,csys,units,varargin)
% ------------------------------------------------------------------------|
% DESCRIPTION:
% [v vx vy vz vxx vxy vxz vyy vyz vzz slay sss] = newton_litho1d0_ver1( ...
%           lon,lat,density,R1,R2,dSp,G,fntl,ngl,lays,nres,Rh,tlon,tlat,...
%           tg,ief,iefc,layerss,csys,units)
% - Integrates a spherical triangular layered model of the density distr.
% ------------------------------------------------------------------------|
% INPUT:
% lon - latitude of density nodes, [deg]
% lat - latitude of density nodes, [deg]
% density - density, [kg/m3]
% R1 - lower boundary wrt to layer, [m]
% R2 - upper boundary wrt to layer, [m]
% dSp - point-wise differential surface - for each point, [sterad]
%     - represents the area of the spherical surface near the point
%     - its sum over the global should tend to 4*pi
% G - gravitational constant, [m3/s2]
% fntl - specifies what quantities to be computed
%      - 'all' - potential, gravitational vector and its cart. gradients
%              - most time consuming!!!
%      - 'v' - potential
%      - 'vi' - gravitational vector
%      - 'vij' - gravitational tensor
% ngl - number of G-L grid points used for the mean kernel computation
%     - NOT WORKING YET, set ngl=1 ... G-L will not be used
%     - in terms of SH degree, e.g. ngl = 10
% lays - [layer1 layer2] computes integral for all layers in layer1:layer2
%      - can be a vector of indices, where indices correspond with the
%        row dimension of [R1,R2,density]. Applies if length(lays)>2
% nres - resolution of the target grid in terms of spherical harmonic degree
% Rh - radius of the sphere with the target grid, [m]
%    - can be a vector for scattered points
% tlon - longitude of the target grid, [tlon1 tlon2] in [deg]
%      - if [-180 180] ... grid will be set to global [nres+1,2*nres]
%      - for a vector with dim>2, the vector will be used
%      - vector with longitude for tg = 'scatter', [deg]
% tlat - latitude of the target grid, [tlat1(north) tlat2(south)] in [deg]
%      - if [90 -90] ... grid will be set to global with nres+1 points
%      - for a vector with dim>2, the vector will be used for scattered
%        data
%      - vector with latitude for tg = 'scatter', [deg]
% tg - mode for the target grid
%    - 'glg' - Gauss-Legendre grid wrt nres
%    - 'eag' - equi-angular grid wrt nres
%    - scatter - scatter point representation
%              - it must hold that dim of lon, lat, Rh is the same
% ief - mode for including/excluding data from the model around one
%       particular point.
%     - 'in' - set to include only data up to a specific radius (erad)
%            - computes the effect of masses out of the target region
%     - 'ex' - set to exclude all data up to a specific radius (erad)
%            - computes the effect of masses out of the target region
% iefs - ief settings, use [elat elon erad], where
%      - elon - longitude of a center point, doesn't have to be from "lon"
%      - elat - latitude of a center point, doesn't have to be from "lat"
%      - erad - radius for including/excluding data points inside 
%               the rectangle [elat,elon] + erad
%             - can be easily changed to psi ... modification needed!!!
%      -set iefs = [0 0 0] not to use this option
% layerss - string cell with names of all layers
% csys - a string defining the output coordinate system 
%         - 'eff'/'EFF' - global coor. system (xy in the equatorial plane)
%         - 'lnof'/'LNOF', 'nwu'/'NWU' - local system (north-west-up)
% units - output units - potential is always left in SI
%       - case insensitive, use:
%           -'si'/'SI'/'Si' ... meter, second
%           -'grav'/'GRAV'/'Grav' ... mGal, E
% ------------------------------------------------------------------------|
% OUTPUT:
% v - gravitational potential, [m2/s2], 3D array for >1 layers
% v* - gravitational vector, 3D array for >1 layers
% v** - gravitational tensor, 3D array for >1 layers
% slay - list of layers
% sss - a string for given settings
% ------------------------------------------------------------------------|
% OTHERS:
% JS(05-2021)
% - needs functions glgrid, eff2lnof*
% ------------------------------------------------------------------------|
% REFERENCES
% ------------------------------------------------------------------------|
% clear all;
% % % sim input
% 
% % % % example of settings
% H = 150e3;
% Rh = 6371e3 + H; % above this sphere
% G = 6.67408e-11;
% 
% % % % target grid
% nres = 180; % resolution in terms of SH degree
% ngl = 1; % Gauss-Legendre resolution for mean kernel computation
% 
% % % % target area ... for global use -90/90, -180/180-180/nres
% tlon = [-30 30];%[-180 180-180/nres];
% tlat = [-30 30];
% tg = 'glg';
% 
% % % % rectangle exclude/include from the integration ... in degree
% ief = 'in'; % in ex
% iefs = [0 0 0];
% 
% lays = [128 154];
% units = 'grav';
% csys = 'lnof';
% fntl = 'all';
%%
% % % ====================================================================
% % % ====================================================================
% % % ====================================================================

units = lower(units);
switch units
    case 'si'
        ui = 1;
        uij = 1;
    case 'grav'
        ui = 1e5;
        uij = 1e9;
end

ro = 180/pi;

% % % % GL points for mean kernel 
% [thmk, accur, wimk, swi] = glgrid(ngl);
% thmk = flipdim(thmk,1); % because it will be needed below
% clear accur swi;

% % % ====================================================================
% % % Coordinates of the target grid
% % % ====================================================================

if length(tg)==3 % global or local grid
    switch tg
        case  'glg'
            [thzeros, accur, wi, swi] = glgrid(nres+1);
            latv0 = 90-thzeros;
            % % % find those within the target area
            ilat = find(latv0>=tlat(2) & latv0<=tlat(1)); %south to north
            latv = latv0(ilat);
        case 'eag'
            latv0 = -(tlat(1):180/nres:tlat(2));
            ilat = find(latv0>=tlat(1) & latv0<=tlat(2));
            latv = latv0(ilat);
    end
    if tlon(1)==-180 && (tlon(2)==180 || tlon(2)==180-180/nres)  % to decide if it is global or local grid
        dlon = 180/nres;
    elseif tlon(1)==0 && tlon(2)==360
        dlon = 0;
    else
        dlon = 0;
    end
    lonv = (tlon(1):180/nres:tlon(2)-dlon);
    
    % % % target grid
    [lont,latt] = meshgrid(lonv,latv);
    [st1,st2] = size(lont);
    
    sRh = num2str(Rh);
    snres=['r' num2str(nres)];
    
else % vector, can contain scatter points
    tg = 'scatter';
    lont = tlon;
    latt = tlat;
   
    if length(lont)~=length(latt) || length(Rh)~=length(latt)
       error('All three vectors should have the same length'); 
    end
    [st1, st2] = size(lont);
    Rh = reshape(Rh,st1*st2,1);
    sRh = 'Rh-user';
    snres = 'r-user';
end

xa = Rh.*reshape(cos(latt/ro).*cos(lont/ro),st1*st2,1);
ya = Rh.*reshape(cos(latt/ro).*sin(lont/ro),st1*st2,1);
za = Rh.*reshape(sin(latt/ro),st1*st2,1);

%%
% % % ====================================================================
% % % constraining input
% % % ====================================================================

elon = iefs(1);
elat = iefs(2);
erad = iefs(3);

ie0 = find(abs(elat-lat(1,:)) < erad & abs(elon-lon(1,:)) < erad); % find all in the rectangle
ie = (1:length(dSp)); % all nodes
if strcmp(ief,'ex') == 1 && erad~=0 % exclude
    ie(ie0) = [];
elseif strcmp(ief,'in') == 1 && erad~=0 % include
    ie = ie0; 
end
if elat==0 && elon==0 && erad==0
   ief = '';    
end
% figure; scatter(lon(1,ie),lat(1,ie),'.');
% dfsdf
% % % layers to be used
if size(lays,1)*size(lays,2)==2
   [mnlay] = min(lays);
   [mxlay] = max(lays);
   ilay = mnlay:mxlay; 
else
   ilay = lays;
end

slay = layerss(ilay);

if length(varargin)>0
   lonp = varargin{1}; 
   latp = varargin{2};
   
    [s1p,s2p] = size(lonp);
   if s2p ~= size(lon,2)
       error('For varargin, size(lonp,2) must be equal size(lon,2)');
   end
   lonp=lonp(:,ie);
   latp=latp(:,ie);
end

% % %  take nodes for particular layers and particular input area
dSp = dSp(1,ie);
density = density(ilay,ie);
lat = lat(:,ie);
lon = lon(:,ie);
R1 = R1(ilay,ie);
R2 = R2(ilay,ie);

[s1,s2] = size(density);

% % % volume of the element for each nodal point
dV = 1/3*repmat(dSp,s1,1).*(R2.^3 - R1.^3); % [sum(dV(:)) % sum(dSp)]
rhodV = density.*dV;

% % % find NaN indices
ii = find(isnan(rhodV)==1);

% % % correct density and R1, R2 because of possible NaNs
rhodV(ii) = 0;
R1(ii) = 6200000; % can be whatever because rhodV is zero there
R2(ii) = 6200000;

% % Matrix with the radius at G-L points
% Rm = zeros(s1,s2,ngl);
% for i=1:ngl % for ngl=1, it results to (R1+R2)/2;
%     Rm(:,:,i) = ((R2-R1)*cos(thmk(i)/ro) + R1+R2)/2;
% end
% % % rewrites above
Rm =  (R2 + R1)/2;

% clear W;
% W(1,1,:) = wimk;
% W = repmat(W,[s1 s2 1]);

% % % integration
tic;
if length(varargin)>0
    
% merge the central point with the surrounding     
    lon=[lon;lonp]/ro;
    lat=[lat;latp]/ro;
    
% neglecting surrounding points
%     lon=repmat(lon,s1p+1,1)/ro;
%     lat=repmat(lat,s1p+1,1)/ro;
       
    % % % input coordinates - input grid    
%     Rmm= permute(repmat(Rm,1,1,s1p+1),[3 2 1]); 
    Rmm= repmat(Rm,1,1,s1p+1); 
    
    lat=permute(repmat(lat',1,1,s1),[3 1 2]);
    lon=permute(repmat(lon',1,1,s1),[3 1 2]);
    
    xL = Rmm.*cos(lat).*cos(lon);
    yL = Rmm.*cos(lat).*sin(lon);
    zL = Rmm.*sin(lat);

    xL(isnan(xL))=0;
    yL(isnan(yL))=0;
    zL(isnan(zL))=0;
    
    fntl = lower(fntl);
    switch fntl

        case 'v'
            w0=[s1p;ones(s1p,1)];
            swi=sum(w0);
            wi = permute(repmat(w0,1,s1,s2),[2 3 1]);
            v = zeros(s1,st1*st2); % % % allocation of output
            for k=1:st1*st2
                K = squeeze(sum(wi./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2),3)/2)/swi;
                v(:,k)=sum(rhodV.*K,2); % integration
            end
            v = G*permute((reshape(v',st1,st2,[])),[1 2 3]);
            vx = 0; vy = 0; vz = 0;
            vxx = 0; vxy = 0; vxz = 0; vyy = 0; vyz = 0; vzz = 0;
            
        case 'vi'
            w0=[6.5*s1p;ones(s1p,1)];
            swi=sum(w0);
            wi = permute(repmat(w0,1,s1,s2),[2 3 1]);
            vx = zeros(s1,st1*st2); % % % allocation of output
            vy = zeros(s1,st1*st2);
            vz = zeros(s1,st1*st2);
                        
            for k=1:st1*st2
                il = 1./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2); 
                %             K = sum(W.*(il.*il.*il.*(xa(k) - xL)),3)/2;  % x  ... divide by 2 if u use GL
                K = sum(wi.*(il.*il.*il.*(xa(k) - xL)),3)/swi;  % x
                vx(:,k) = sum(rhodV.*K,2);
                K = sum(wi.*(il.*il.*il.*(ya(k) - yL)),3)/swi;  % y
                vy(:,k) = sum(rhodV.*K,2);
                K = sum(wi.*(il.*il.*il.*(za(k) - zL)),3)/swi;  % z
                vz(:,k) = sum(rhodV.*K,2);
            end
            vx = ui*G*permute((reshape(vx',st1,st2,[])),[1 2 3]); % transpose 3D matrix
            vy = ui*G*permute((reshape(vy',st1,st2,[])),[1 2 3]);
            vz = ui*G*permute((reshape(vz',st1,st2,[])),[1 2 3]);
            v = 0;
            vxx = 0; vxy = 0; vxz = 0; vyy = 0; vyz = 0; vzz = 0;
            
    end
    
else % single-point kernel
    
    % % % unit coordinates - input grid
    xL = Rm.*repmat(cos(lat(1,:)/ro).*cos(lon(1,:)/ro),s1,1,ngl);
    yL = Rm.*repmat(cos(lat(1,:)/ro).*sin(lon(1,:)/ro),s1,1,ngl);
    zL = Rm.*repmat(sin(lat(1,:)/ro),s1,1,ngl);
    
    fntl = lower(fntl);
    switch fntl
        case 'v'
            v = zeros(s1,st1*st2); % % % allocation of output
            
            for k=1:st1*st2
                %             K = sum(W./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2),3)/2; 
                K = sum(1./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2),3)/2; 
                v(:,k)=sum(rhodV.*K,2); % integration
            end
            
            v = G*permute((reshape(v',st1,st2,[])),[1 2 3]);
            vx = 0; vy = 0; vz = 0;
            vxx = 0; vxy = 0; vxz = 0; vyy = 0; vyz = 0; vzz = 0;
            
        case 'vi'
            vx = zeros(s1,st1*st2); % % % allocation of output
            vy = zeros(s1,st1*st2);
            vz = zeros(s1,st1*st2);
            
            for k=1:st1*st2
                
                il = 1./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2); 
                %             K = sum(W.*(il.*il.*il.*(xa(k) - xL)),3)/2;  % x  ... divide by 2 if u use GL
                K = sum(1.*(il.*il.*il.*(xa(k) - xL)),3);  % x
                vx(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(il.*il.*il.*(ya(k) - yL)),3);  % y
                vy(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(il.*il.*il.*(za(k) - zL)),3);  % z
                vz(:,k) = sum(rhodV.*K,2);
            end
            
            vx = ui*G*permute((reshape(vx',st1,st2,[])),[1 2 3]); % transpose 3D matrix
            vy = ui*G*permute((reshape(vy',st1,st2,[])),[1 2 3]);
            vz = ui*G*permute((reshape(vz',st1,st2,[])),[1 2 3]);
            v = 0;
            vxx = 0; vxy = 0; vxz = 0; vyy = 0; vyz = 0; vzz = 0;
            
        case 'vij'
            vxy = zeros(s1,st1*st2);
            vyz = zeros(s1,st1*st2);
            vxz = zeros(s1,st1*st2);
            vxx = zeros(s1,st1*st2);
            vyy = zeros(s1,st1*st2);
            vzz = zeros(s1,st1*st2);
            
            tic;
            for k=1:st1*st2
                
                il = 1./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2);
                
                %             K = sum(W.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).^2) - il.*il.*il),3)/2;  % ... xx
                K = sum(1.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).^2) - il.*il.*il),3);  % ... xx
                vxx(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).*(ya(k) - yL))),3);    % ... xy
                vxy(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).*(za(k) - zL))),3);    % ... xz
                vxz(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*(ya(k) - yL).^2 - il.*il.*il),3);    % ... yy
                vyy(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*((ya(k) - yL).*(za(k) - zL))),3);    % ... yz
                vyz(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*(za(k) - zL).^2 - il.*il.*il),3);    % ... zz
                vzz(:,k) = sum(rhodV.*K,2);
            end
            
            vxx = uij*G*permute((reshape(vxx',st1,st2,[])),[1 2 3]);
            vxy = uij*G*permute((reshape(vxy',st1,st2,[])),[1 2 3]);
            vxz = uij*G*permute((reshape(vxz',st1,st2,[])),[1 2 3]);
            vyy = uij*G*permute((reshape(vyy',st1,st2,[])),[1 2 3]);
            vyz = uij*G*permute((reshape(vyz',st1,st2,[])),[1 2 3]);
            vzz = uij*G*permute((reshape(vzz',st1,st2,[])),[1 2 3]);
            v = 0;
            vx = 0; vy = 0; vz = 0;
            
        case 'all'
            v = zeros(s1,st1*st2); % % % allocation of output
            vx = zeros(s1,st1*st2);
            vy = zeros(s1,st1*st2);
            vz = zeros(s1,st1*st2);
            vxy = zeros(s1,st1*st2);
            vyz = zeros(s1,st1*st2);
            vxz = zeros(s1,st1*st2);
            vxx = zeros(s1,st1*st2);
            vyy = zeros(s1,st1*st2);
            vzz = zeros(s1,st1*st2);
            
            tic;
            for k=1:st1*st2
                
                il = 1./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2);
                
                %             K = sum(W./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2),3)/2; 
                K = sum(1./sqrt((xa(k) - xL).^2 + (ya(k) - yL).^2 + (za(k) - zL).^2),3);
                v(:,k)=sum(rhodV.*K,2); % integration
                
                K = sum(1.*(il.*il.*il.*(xa(k) - xL)),3);  % x
                vx(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(il.*il.*il.*(ya(k) - yL)),3);  % y
                vy(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(il.*il.*il.*(za(k) - zL)),3);  % z
                vz(:,k) = sum(rhodV.*K,2);
                
                K = sum(1.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).^2) - il.*il.*il),3);  % ... xx
                vxx(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).*(ya(k) - yL))),3);    % ... xy
                vxy(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*((xa(k) - xL).*(za(k) - zL))),3);    % ... xz
                vxz(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*(ya(k) - yL).^2 - il.*il.*il),3);    % ... yy
                vyy(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*((ya(k) - yL).*(za(k) - zL))),3);    % ... yz
                vyz(:,k) = sum(rhodV.*K,2);
                K = sum(1.*(3*il.*il.*il.*il.*il.*(za(k) - zL).^2 - il.*il.*il),3);    % ... zz
                vzz(:,k) = sum(rhodV.*K,2);
            end
            
            v = ui*G*permute((reshape(v',st1,st2,[])),[1 2 3]);
            
            vx = ui*G*permute((reshape(vx',st1,st2,[])),[1 2 3]);
            vy = ui*G*permute((reshape(vy',st1,st2,[])),[1 2 3]);
            vz = ui*G*permute((reshape(vz',st1,st2,[])),[1 2 3]);
            
            vxx = uij*G*permute((reshape(vxx',st1,st2,[])),[1 2 3]);
            vxy = uij*G*permute((reshape(vxy',st1,st2,[])),[1 2 3]);
            vxz = uij*G*permute((reshape(vxz',st1,st2,[])),[1 2 3]);
            vyy = uij*G*permute((reshape(vyy',st1,st2,[])),[1 2 3]);
            vyz = uij*G*permute((reshape(vyz',st1,st2,[])),[1 2 3]);
            vzz = uij*G*permute((reshape(vzz',st1,st2,[])),[1 2 3]);
    end
end
tt=toc;

% % % rotation into the target frame
if strcmpi(csys,'eff')~=1 && length(tg)==3
    if strcmp(fntl,'vi')==1
        w = zeros(st1,st2,3);
        for i=1:s1
            w(:,:,1) = vx(:,:,i);
            w(:,:,2) = vy(:,:,i);
            w(:,:,3) = vz(:,:,i);
            [w] = eff2lnofe_vi(w,(90-latt)/ro,lont/ro,latt/ro);
            vx(:,:,i) = w(:,:,1);
            vy(:,:,i) = w(:,:,2);
            vz(:,:,i) = w(:,:,3);
        end
    elseif strcmp(fntl,'vij')==1
        ww = zeros(st1,st2,6);
        for i=1:s1
            ww(:,:,1) = vxx(:,:,i);
            ww(:,:,2) = vxy(:,:,i);
            ww(:,:,3) = vxz(:,:,i);
            ww(:,:,4) = vyy(:,:,i);
            ww(:,:,5) = vyz(:,:,i);
            ww(:,:,6) = vzz(:,:,i);            
            [ww] = eff2lnofe_vij(ww,(90-latt)/ro,lont/ro,latt/ro);
            vxx(:,:,i) = ww(:,:,1);
            vxy(:,:,i) = ww(:,:,2);
            vxz(:,:,i) = ww(:,:,3);
            vyy(:,:,i) = ww(:,:,4);
            vyz(:,:,i) = ww(:,:,5);
            vzz(:,:,i) = ww(:,:,6);
        end
        
    elseif  strcmpi(fntl,'all')==1
        w = zeros(st1,st2,3);
        ww= zeros(st1,st2,6);
        for i=1:s1   
            w(:,:,1) = vx(:,:,i);
            w(:,:,2) = vy(:,:,i);
            w(:,:,3) = vz(:,:,i);
            ww(:,:,1) = vxx(:,:,i);
            ww(:,:,2) = vxy(:,:,i);
            ww(:,:,3) = vxz(:,:,i);
            ww(:,:,4) = vyy(:,:,i);
            ww(:,:,5) = vyz(:,:,i);
            ww(:,:,6) = vzz(:,:,i);     
         
            [w,ww] = eff2lnofe(w,ww,pi/2-latt/ro,lont/ro,latt/ro);
            vx(:,:,i) = w(:,:,1);
            vy(:,:,i) = w(:,:,2);
            vz(:,:,i) = w(:,:,3);
            vxx(:,:,i) = ww(:,:,1);
            vxy(:,:,i) = ww(:,:,2);
            vxz(:,:,i) = ww(:,:,3);
            vyy(:,:,i) = ww(:,:,4);
            vyz(:,:,i) = ww(:,:,5);
            vzz(:,:,i) = ww(:,:,6);
        end     
    end
end

if length(lays)==1
    layer1 = lays(1);
    layer2 = lays(1);
    ssss=[num2str(layer1) '-' num2str(layer2)];
elseif length(lays)==2
    layer1 = lays(1);
    layer2 = lays(2);
    ssss=[num2str(layer1) '-' num2str(layer2)];
elseif length(lays)>=2 && sum(abs(diff(lays))) == length(lays)-1 % layer by layer
    layer1 = lays(1);
    layer2 = lays(1);
    ssss=[num2str(layer1) '-' num2str(layer2)];
else length(lays)>=2 && sum(abs(diff(lays))) ~= length(lays)-1 % random layers
    ssss=[num2str(layer1) 'Frand-' num2str(layer2) 'Lrand']; % First and last
end

sss = [sRh '_' snres '_'  csys '_' fntl '_units-' units '_grid-' tg '_layers' ssss '_' ief num2str(elon) '-' num2str(elat) '-' num2str(erad)];

