function m = upgrademesh( m, upgradeto )
%m = upgrademesh( m )
%   Ensure the mesh M is compatible with the current version of the software.

    if ~isfield( m, 'versioninfo' )
        m = upgradeVeryOldMesh( m );
        % m now has version number 0.
    end
    
    setGlobals();
    global gMAX_MESHVERSION gOLD_NUMRESERVEDMGENS gNEW_NUMRESERVEDMGENS
    global gOLD_K_POL
    global gNEW_POLARISER
    global gOBSOLETEFIELDS gOBSOLETESTATICPROPS gDEFAULTFIELDS
    global gOurViewParams

    if isfield( m.versioninfo, 'generalversion' )
        m.versioninfo = struct( 'meshversion', 0, ...
            'mgenversion', m.versioninfo.generalversion );
    end
    if m.versioninfo.mgenversion==0
        pol_mgen = gOLD_K_POL;
    else
        pol_mgen = gNEW_POLARISER;
    end
    global gGlobalProps gGlobalDynamicProps
    global gDefaultPlotOptions;
    
    if nargin < 2
        upgradeto = gMAX_MESHVERSION;
    end
    
    numMorphogens = size(m.morphogens,2);

        
    % 2008 Jan 11  Default model name is now 'untitled' instead of the empty string.
    if strcmp( m.globalProps.modelname, 'untitled' )
        m.globalProps.modelname = '';
    end
    
    m.globalProps = replacefields( m.globalProps, ...
        'bioAsplitproc', 'bioApresplitproc', ...
        'bioAsplatproc', 'bioApostsplitproc' );

    if isfield( m, 'canvascolor' )
        m.plotdefaults.canvascolor = m.canvascolor;
        m = rmfield( m, 'canvascolor' );
    end
    
    m.plotdefaults = upgradePlotoptions( m.plotdefaults );

    if isfield( m.globalProps, 'defaultViewParams' )
        m.globalProps.defaultViewParams = replacefields( ...
            m.globalProps.defaultViewParams, 'CameraUp', 'CameraUpVector' );
    else
        m.globalProps.defaultViewParams = m.plotdefaults.matlabViewParams;
    end

    if (~isfield( m, 'mgen_transportable' )) || (length(m.mgen_transportable) ~= numMorphogens)
        m.mgen_transportable = false( 1, numMorphogens );
    end

    if ~isfield( m, 'waypoints' )
        m.waypoints = [];
    end

    if ~isfield( m, 'moviescripts' )
        m.moviescripts = [];
    end

    if ~isfield( m, 'movieselected' )
        m.movieselected = 0;
    end
    
    if ~isfield( m, 'transportfield' )
        m.transportfield = cell( 1, numMorphogens );
    end

    if ~isfield( m, 'cellFrames' )
        if isfield( m, 'celldata') && isfield( m.celldata, 'cellFrame')
            m.cellFrames = reshape( [m.celldata.cellFrame], 3, 3, [] );
            m.celldata = rmfield( m.celldata, 'cellFrame' );
        else
            m.cellFrames = [];
        end
    end

    if ~isfield( m, 'outputcolors' )
        m.outputcolors = struct();
    end
    ocfns = fieldnames( m.outputcolors );
    for i=1:length(ocfns)
        if ~isempty( regexp( ocfns{i}, '^actual' , 'once' ) )
            newocfn = regexprep( ocfns{i}, '^actual', 'resultant' );
            m.outputcolors.(newocfn) = m.outputcolors.(ocfns{i});
            m.outputcolors = rmfield( m.outputcolors, ocfns{i} );
        end
    end
    m.outputcolors = defaultFromStruct( m.outputcolors, gDEFAULTFIELDS.outputcolors );
    

    if isfield( m.globalProps, 'dointernalrotation' )
        m.globalProps = rmfield( m.globalProps, 'dointernalrotation' );
        m.globalProps.performinternalrotation = false;
    end
    if isfield( m.globalProps, 'physicalThickness' )
        if m.globalProps.physicalThickness
            m.globalProps.thicknessMode = 'physical';
        else
            m.globalProps.thicknessMode = 'scaled';
        end
        m.globalProps = rmfield( m.globalProps, 'physicalThickness' );
    end
    if strcmp( m.globalProps.thicknessMode, 'anticurl' )
        m.globalProps.thicknessMode = 'direct';
    end
    if isfield( m.globalProps, 'legend' )
        m.globalProps.legendTemplate = gGlobalProps.legendTemplate;
        if ~isempty( m.globalProps.legend )
            m.globalProps.legendTemplate = [ m.globalProps.legendTemplate ...
                                             ' ' ...
                                             m.globalProps.legend ];
        end
        m.globalProps = rmfield( m.globalProps, 'legend' );
    end
    
    if isfield( m, 'invisible' )
        if ~isfield( m, 'visible' )
            m = findVisiblePart( m );
        end
        m = rmfield( m, 'invisible' );
    end

    % Remove obsolete fields.
    m = safermfield( m, gOBSOLETEFIELDS );
    m.globalProps = safermfield( m.globalProps, gOBSOLETESTATICPROPS );
    
    m = defaultFromStruct( m, gDEFAULTFIELDS );
    
    if iscell( m.stagetimes )
        m.stagetimes = [];
    end
    
    m.meshparams = defaultFromStruct( m.meshparams, ...
                                      struct( ...
                                        'type', 'unknown', ...
                                        'randomness', 0 ) );
    
    if isfield( m, 'mgencolors' )
        m.mgenposcolors = m.mgencolors;
        m = rmfield( m, 'mgencolors' );
    end
    if ~isfield( m, 'mgenposcolors' )
        m.mgenposcolors = HSVtoRGB( [ (0:1:(numMorphogens-1))'/12, ones( numMorphogens, 2 ) ] )';
    end
    if ~isfield( m, 'mgennegcolors' )
        m.mgennegcolors = oppositeColor( m.mgenposcolors' )';
    end
    
    if ~isfield( m, 'mgen_plotpriority' )
        m.mgen_plotpriority = zeros( 1, numMorphogens );
    end

    if ~isfield( m, 'mgen_plotthreshold' )
        m.mgen_plotthreshold = zeros( 1, numMorphogens );
    end
    
    if isfield( m.secondlayer, 'cells' ) && (~isfield( m.secondlayer, 'side' ))
        m.secondlayer = true( length( m.secondlayer.cells ) );
    end

    if all(m.displacements==0)
        m.displacements = [];
    end
    
    gdpfn = fieldnames( gGlobalDynamicProps );
    if ~isfield( m, 'globalDynamicProps' )
        gdp = gGlobalDynamicProps;
        for i=1:length(gdpfn)
            fn = gdpfn{i};
            if isfield( m.globalProps, fn )
                gdp.(fn) = m.globalProps.(fn);
            end
        end
        m.globalProps = safermfield( m.globalProps, gdpfn );
        m.globalDynamicProps = gdp;
    else
        for i=1:length(gdpfn)
            if isfield( m.globalProps, gdpfn{i} )
                m.globalDynamicProps.(gdpfn{i}) = m.globalProps.(gdpfn{i});
                m.globalProps = rmfield( m.globalProps, gdpfn{i} );
            end
        end
    end
    if ~isfield( m.globalDynamicProps, 'staticreadonly' )
        m.globalDynamicProps.staticreadonly = false;
    end

    m.globalProps.flatten = false;  % Obsolete, to be deleted in a future version.
    m.globalProps.flattenratio = 1; % Obsolete, to be deleted in a future version.
    if ~isfield( m.globalProps, 'relativepolgrad' )
         m.globalProps.relativepolgrad = gGlobalProps.relativepolgrad;
    end
    if ~isfield( m.globalProps, 'mingradient' )
         m.globalProps.mingradient = gGlobalProps.mingradient;
    end
    if ~isfield( m.globalProps, 'userpolarisation' )
         m.globalProps.userpolarisation = gGlobalProps.userpolarisation;
    end
    if ~isfield( m.globalProps, 'usepolfreezebc' )
         m.globalProps.usepolfreezebc = gGlobalProps.usepolfreezebc;
    elseif ischar( m.globalProps.usepolfreezebc )
        m.globalProps.usepolfreezebc = gGlobalProps.usepolfreezebc;
    end
    if ~isfield( m.globalProps, 'usefrozengradient' )
         m.globalProps.usefrozengradient = gGlobalProps.usefrozengradient;
    end
    
    m = inflatemesh( m );
    
    % Use integers, where integers are what is wanted.
    % The field-testing here is because we want to be able to handle stripped
    % meshes as well.
    m.tricellvxs = forcetype( m.tricellvxs, 'int32' );
    if isfield( m, 'edgecells' )
        m.edgecells = forcetype( m.edgecells, 'int32' );
    end
    if isfield( m, 'celledges' )
        m.celledges = forcetype( m.celledges, 'int32' );
    end
    if hasSecondLayer( m )
        if isfield( m.secondlayer, 'edges' )
            m.secondlayer.edges = forcetype( m.secondlayer.edges, 'int32' );
        end
        m.secondlayer.cloneindex = forcetype( m.secondlayer.cloneindex, 'int32' );
        if isfield( m.secondlayer.cells, 'edges' )
            for i=1:length(m.secondlayer.cells)
                m.secondlayer.cells(i).edges = forcetype( m.secondlayer.cells(i).edges, 'int32' );
            end
        end
    end
    if isfield( m, 'nodecelledges' )
        for i=1:length(m.nodecelledges)
            m.nodecelledges{i} = int32( m.nodecelledges{i} );
        end
    end
    fns = fieldnames( m.mgenNameToIndex );
    for i=1:length(fns)
        m.mgenNameToIndex.(fns{i}) = int32( m.mgenNameToIndex.(fns{i}) );
    end
    
    
    if ~isfield( m.globalProps, 'starttime' )
        m.globalProps.starttime = ...
            m.globalDynamicProps.currenttime ...
            - m.globalProps.timestep * m.globalDynamicProps.currentIter;
    end
    if ~isfield( m.globalProps, 'unitbulkmodulus' )
        if m.globalProps.bulkmodulus ~= 1
            m = leaf_setproperty( m, 'bulkmodulus', 1 );
        end
    end
    m.globalProps = defaultFromStruct( m.globalProps, gGlobalProps );
    m.globalProps = safermfield( m.globalProps, 'MESHVERSION' );
    if isnumeric(m.globalProps.defaultinterp)
        switch m.globalProps.defaultinterp
            case 1
                m.globalProps.defaultinterp = 'mid';
            case 2
                m.globalProps.defaultinterp = 'min';
            case 3
                m.globalProps.defaultinterp = 'max';
            otherwise
        end
    end
    m.globalProps.coderevision = int32(m.globalProps.coderevision);
    m.globalProps.maxFEcells = int32(m.globalProps.maxFEcells);
    m.globalProps.inittotalcells = int32(m.globalProps.inittotalcells);
    m.globalProps.maxBioAcells = int32(m.globalProps.maxBioAcells);
    m.globalProps.canceldrift = logical(m.globalProps.canceldrift);
    m.globalProps.allowInteraction = logical(m.globalProps.allowInteraction);
    if ~isfield( m, 'mgen_interpType' )
        if m.versioninfo.mgenversion==0
            numStdMgens = gOLD_NUMRESERVEDMGENS;
        else
            numStdMgens = gNEW_NUMRESERVEDMGENS;
        end
        for i=1:numStdMgens
            m.mgen_interpType{i} = 'mid';
        end
        for i=numStdMgens+1:numMorphogens
            m.mgen_interpType{i} = m.globalProps.defaultinterp;
        end
    elseif ~iscell( m.mgen_interpType )
        numMorphogens = length(m.mgen_interpType);
        newInterpType = cell(1,numMorphogens);
        for i=1:numMorphogens
            switch m.mgen_interpType(i)
                case 2
                    newInterpType{i} = 'min';
                case 3
                    newInterpType{i} = 'max';
                otherwise
                    newInterpType{i} = 'mid';
            end
        end
        m.mgen_interpType = newInterpType;
    end
    
    if ~isfield( m, 'mgenswitch' )
        m.mgenswitch = ones( 1, size(m.mutantLevel,2) );
    end
    
    if isfield( m.globalProps, 'laststage' )
        m.globalDynamicProps.laststagesuffix = '';
        m.globalProps = rmfield( m.globalProps, 'laststage' );
    end
    
    m = updateElasticity( m );
    
    if hasSecondLayer( m )
        if isfield( m.secondlayer, 'colors' )
            m = setSecondLayerColorInfo( m, m.secondlayer.colors, m.secondlayer.colorvariation );
            m.secondlayer = rmfield( m.secondlayer, { 'colors', 'colorvariation', 'colorparams' } );
        end
    else
        m = setSecondLayerColorInfo( m );
    end
    if size(m.globalProps.colors,1)==1
        m.globalProps.colors = [ m.globalProps.colors; 1-m.globalProps.colors ];
    end

    if ~isfield( m, 'mgen_production' )
        m.mgen_production = zeros(size(m.morphogens));
    end
    
    % Remove obsolete fields of m.celldata.
    m.celldata = safermfield( m.celldata, ...
        'cellThermDiffLocalTensor', ...
        'cellThermDiffDir', ...
        'cellThermDiffGlobalTensor', ...
        'cellThermExpLocalTensor', ...
        'cellThermExpDir', ...
        'majorThermExp', ...
        'minorThermExp', ...
        'avResidStrain', ...
        'residualStress', ...
        'normResidStrain', ...
        'normResidStress', ...
        'actualGrowthParams' );
    if isfield(m.celldata, 'residualStrain') && (size(m.celldata(1).residualStrain,2)==1)
        for i=1:length(m.celldata)
            m.celldata(i).residualStrain = repmat( m.celldata(i).residualStrain, 1, 6 );
        end
    end
    if isfield(m.celldata, 'residualStress') && (size(m.celldata(1).residualStress,2)==1)
        for i=1:length(m.celldata)
            m.celldata(i).residualStress = repmat( m.celldata(i).residualStress, 1, 6 );
        end
    end
    
    % Older versions of GFtbox could leave m.outputs in an inconsistent
    % state (e.g. if leaf_load were used to replace the mesh by another).
    % Here we test to see if every component of m.outputs is the size it
    % should be, and if not, discard and recompute it.
    numcells = size( m.tricellvxs, 1 );
    goodOutputs = isfield( m, 'outputs' ) ...
        && (size( m.outputs.specifiedstrain.A, 1 ) == numcells) ...
        && (size( m.outputs.specifiedstrain.B, 1 ) == numcells) ...
        && (size( m.outputs.actualstrain.A, 1 ) == numcells) ...
        && (size( m.outputs.actualstrain.B, 1 ) == numcells) ...
        && (size( m.outputs.residualstrain.A, 1 ) == numcells) ...
        && (size( m.outputs.residualstrain.B, 1 ) == numcells) ...
        && (size( m.outputs.rotations, 3 ) == numcells);
    if ~goodOutputs
        m = calculateOutputs( m );
    end
    
    if isfield( m, 'cellareas' ) && (size(m.cellareas,2) > 1)
        m.cellareas = m.cellareas(:);
    end
    
    if ~isfield( m, 'rewriteIFneeded' )
        m.rewriteIFneeded = false;
    end

    PERCELL_ELASTICITY_IS_UNIMPLEMENTED = true;
    if PERCELL_ELASTICITY_IS_UNIMPLEMENTED || ~isfield( m, 'cellstiffness' )
        m.cellbulkmodulus = repmat( m.globalProps.bulkmodulus, [size(m.tricellvxs,1), 1] );
        m.cellpoisson = repmat( m.globalProps.poissonsRatio, [size(m.tricellvxs,1), 1] );
        m.cellstiffness = repmat( m.globalProps.D, [1,1,size(m.tricellvxs,1)] );
    end

    if ~isfield( m, 'polfreeze' )
        m.polfreeze = reshape( m.morphogens( m.tricellvxs, pol_mgen ), size(m.tricellvxs) );
    end
    
    if ~isfield( m, 'polfreezebc' )
        for ci=1:size(m.polfreeze,1)
            m.polfreezebc(ci,:) = vec2bc( m.gradpolgrowth(ci,:), ...
                                          m.nodes( m.tricellvxs(ci,:), : ) );
        end
    end
    
    if ~isfield( m, 'polfrozen' )
        m.polfrozen = false( size(m.tricellvxs,1), 1 );
    elseif size( m.polfrozen, 1 )==1
        m.polfrozen = m.polfrozen(:);
    end
    
  % if ~isfield( m, 'meshparams' )
  %     m.meshparams = struct( 'type', 'unknown' );
  % end
    
    m.globalDynamicProps = defaultFromStruct( m.globalDynamicProps, gGlobalDynamicProps  );
    
% 2007 Jun 12
    [ok,m] = validmesh( m );
    if ~ok
        % Nothing -- validmesh will already have written error messages.
    end

    if m.versioninfo.meshversion >= upgradeto
        return;
    end
end

function x = forcetype( x, fieldtype )
    if ~strcmp( class( x ), fieldtype )
        fieldtypehandle = str2func(fieldtype);
        x = fieldtypehandle(x);
    end
end

