function m = rescaleSpaceTime( m, spaceunitname, spacescale, timeunitname, timescale )
%m = rescaleSpaceTime( m, spaceunitname, spacescale, timeunitname, timescale )
%   Choose new units for space and time.  m is physically unchanged by
%   this, and therefore its growth behaviour should be unaltered.  We could
%   refrain from changing m at all, and merely store the factors that
%   convert between internal units and units displayed to the user.
%   However, it is useful to actually rescale m itself, as a check on the
%   invariance of our calculations to such rescaling.
%
%   The new space unit is spacescale old units.  Therefore all of the
%   numbers that represent distances or positions must be divided by
%   spacescale.  Areas must be divided by spacescale^2.  Gradients must be
%   multiplied by spacescale.
%
%   The new time unit is timescale old units.  Therefore every
%   numerical time must be divided by timescale.  Every rate of change must
%   be multiplied by timescale.  Every decay rate must be raised to the
%   power timescale.  Every growth rate g should in principle be replaced
%   by (1+g)^timescale - 1; when g is small, this is close to g*timescale.
%
%   Elasticity coefficients are scale-free.
%
%   Diffusion constants must be multiplied by timescale/spacescale^2.

    setGlobals();
    global gOLD_STRAINRET;
    global gNEW_STRAINRET;
    if m.versioninfo.mgenversion==0
        STRAINRET_MGEN = gOLD_STRAINRET;
    else
        STRAINRET_MGEN = gNEW_STRAINRET;
    end
    
    if ~isempty(spaceunitname)
        m.globalProps.distunitname = spaceunitname;
    end
    if ~isempty(timeunitname)
        m.globalProps.timeunitname = timeunitname;
    end

    if spacescale ~= 1
        spacescalesq = spacescale*spacescale;
        m.nodes = m.nodes/spacescale;
        m.prismnodes = m.prismnodes/spacescale;
        m.displacements = m.displacements/spacescale;
        m.gradpolgrowth = m.gradpolgrowth * spacescale;
        m.cellareas = m.cellareas / spacescalesq;
        m.globalProps.thresholdsq = m.globalProps.thresholdsq / spacescalesq;
        m.globalDynamicProps.cellscale = m.globalDynamicProps.cellscale / spacescale;
        m.globalDynamicProps.thicknessAbsolute = m.globalDynamicProps.thicknessAbsolute / spacescale;
        m.globalDynamicProps.currentArea = m.globalDynamicProps.currentArea / spacescalesq;
        m.globalDynamicProps.previousArea = m.globalDynamicProps.previousArea / spacescalesq;
        m.globalProps.initialArea = m.globalProps.initialArea / spacescalesq;
        m.globalProps.bendunitlength = m.globalProps.bendunitlength / spacescalesq;
        m.plotdefaults.axisRange = m.plotdefaults.axisRange / spacescale;
        m.plotdefaults.clippingDistance = m.plotdefaults.clippingDistance / spacescale;
        m = calcCloneVxCoords( m );
        if hasSecondLayer( m )
            m.secondlayer.splitThreshold = m.secondlayer.splitThreshold/spacescale;
            m.secondlayer.cellarea = m.secondlayer.cellarea/spacescalesq;
            m.secondlayer.celltargetarea = m.secondlayer.celltargetarea/spacescalesq;
            m.secondlayer.averagetargetarea = m.secondlayer.averagetargetarea/spacescalesq;
        end
    end
    
    if timescale ~= 1
        m.timestep = m.timestep/timescale;
        m.currenttime = m.currenttime/timescale;
        m.mgen_absorption = m.mgen_absorption * timescale;
        m.morphogens(:,STRAINRET_MGEN) = ...
            m.morphogens(:,STRAINRET_MGEN) .^ timescale;
    end
    
    Kscale = timescale/(spacescale*spacescale);
    if Kscale ~= 1
        m.absKvector = m.absKvector * Kscale;
    end
end
