function [nc]=nc_add_newcoordsys(nc,varname,RSPstring)

% Assumed that input variable is defined with respect to RSP
% Assumed that the field `nc(:).Attribute(4).Value` (=comment) contains the string `w.r.t. RSP`

%%

if ~exist('RSPstring','var');
    RSPstring='w.r.t. RSP';
    %     RSPstring='t.o.v. RSP';
end

lstr=length(RSPstring);

%%

ifld=find(strcmpi(varname,{nc.Name}));
ifld_orig=ifld;

if isempty(ifld);return;end

if ifld<length(nc)
    nc(ifld+5:end+4)=nc(ifld+1:end);
end
nc(ifld+1:ifld+4)=nc(ifld);


info={...
    'xRD','in RD-coordinates [x]';...
    'yRD','in RD-coordinates [y]';...
    'lat','in WGS84-coordinates [lat]';...
    'lon','in WGS84-coordinates [lon]'...
    };

for i=1:size(info,1);

   ifld = ifld + 1;
   
   nc(ifld).Name                = [nc(ifld_orig).Name '_' info{i,1}];
   nc(ifld).Attribute(1).Value  = [nc(ifld_orig).Attribute(1).Value ' [' info{i,1} ']'];
   
   newstring             = info{i,2};
   lnewstring            = length(newstring);
   comment               = strtrim(nc(ifld_orig).Attribute(4).Value);
   lcomment              = length(comment);
   
   istr                  = strfind(comment,RSPstring);
   lcommentend           = lcomment-(istr+lstr-1);
   if isempty(istr)
       nc(ifld).Attribute(4).Value  = [comment ' ' newstring];
   else
       nc(ifld).Attribute(4).Value = comment;
       nc(ifld).Attribute(4).Value(istr:istr+lnewstring-1) = newstring;
       nc(ifld).Attribute(4).Value(istr+lnewstring:istr+lnewstring+lcommentend-1) = comment(istr+lstr:end);
   end
   
end
