if viz

    % open the file that contains the materials properties:
    fid = fopen(fullfile(pwd,'SWMS_3D.IN','SELECTOR.IN'));
    % skip through to the block containing the materials information:
    fid = goToBlockStart(fid,'B');
    % skip over the block's header:
    fid = skipLines(fid,2);
    % record the position:
    pos = ftell(fid);
    % read how many materials are listed further down:
    C = textscan(fid,'%d',1);
    % assign the number of materials to the variable nMaterials:
    nMaterials = C{1};
    % return to the position previously recorded:
    seekStatus = fseek(fid,pos,'bof');
    % skip two line from the position we just returned to:
    fid = skipLines(fid,2);
    % read the materials properties nMaterials times
    C = textscan(fid,repmat('%f',[1,9]),nMaterials);
    % close the file: 
    closeStatus = fclose(fid);
    % assign each property to its own variable:
    thetaRes.value = C{1};
    thetaRes.units = ['[',swmsLengthUnits,'^{3}_{water}\cdot{}',swmsLengthUnits,'^{-3}_{bulk}]'];
    thetaSat.value = C{2};
    thetaSat.units = ['[',swmsLengthUnits,'^{3}_{water}\cdot{}',swmsLengthUnits,'^{-3}_{bulk}]'];
    parAlpha.value = C{5};
    parAlpha.units = ['[',swmsLengthUnits,'^{-1}]'];
    parN.value = C{6};
    parN.units = ['[-]'];
    clear C

    if strcmp(swmsLengthUnits,'m')
        convFactor = 0.01;
    end



    m = 1-1./parN.value;

    % % % % % % % % % % % % % % % % % % %
    hVec1 = -logspace(-2,3,1000); % cm


    for iMaterial = 1:nMaterials
        A(iMaterial,:) = (1+((parAlpha.value(iMaterial,1)*convFactor.*abs(hVec1)).^parN.value(iMaterial,1))).^m(iMaterial,1);
        thetaAct(iMaterial,:) = thetaRes.value(iMaterial,1)+(thetaSat.value(iMaterial,1)-thetaRes.value(iMaterial,1))./A(iMaterial,:);

        % % % % % % % % % % % % % % % % % % %

        B(iMaterial,:) = (parAlpha.value(iMaterial,:)*convFactor.*abs(hVec1).^(parN.value(iMaterial,:)-1));
        C(iMaterial,:) = (1+((parAlpha.value(iMaterial,:)*convFactor.*abs(hVec1)).^parN.value(iMaterial,:))).^-m(iMaterial,:);
        D(iMaterial,:) = (1+((parAlpha.value(iMaterial,:)*convFactor.*abs(hVec1)).^parN.value(iMaterial,:))).^(m(iMaterial,:)/2);
        relCond(iMaterial,:) = ((1-B(iMaterial,:).*C(iMaterial,:)).^2)./D(iMaterial,:);

        if iMaterial==1
            legStrs{iMaterial} = ['material ',num2str(iMaterial),' (shallow)'];
        elseif iMaterial == nMaterials
            legStrs{iMaterial} = ['material ',num2str(iMaterial),' (deep)'];
        else
            legStrs{iMaterial} = ['material ',num2str(iMaterial)];
        end
    end

    % % % % % % % % % % % % % % % % % % %

    hFig = subplotScreen(1,1,1);
    set(hFig,'numbertitle','off','name','pF',...
        'paperpositionmode','auto','inverthardcopy','off')

    subplot(1,2,1,'fontsize',16)
    plot(thetaAct,hVec1)
    grid on
    set(gca,'yscale','log','xlim',[0,1.0],'ylim',[min(hVec1),max(hVec1)],'ydir','reverse')
    xlabel('water content [m^3/m^3]')
    ylabel('pressure head [cm]')
    legend(legStrs{:})

    subplot(1,2,2,'fontsize',16)
    plot(relCond,hVec1)
    grid on
    set(gca,'xscale','lin','yscale','log','ydir','reverse',...
        'ylim',[min(hVec1),max(hVec1)])
    xlabel('relative conductivity')
    ylabel('pressure head [cm]')
    legend(legStrs{:})


    set(hFig,'paperunits','centimeter')
    pos = get(hFig,'position');
    r = pos(3)/pos(4);
    maxWidth = 26;
    set(hFig,'paperposition',[1,1,maxWidth,maxWidth/r])



    if saveFigs
        print(hFig,'materials.eps','-depsc','-r300','-loose')
        print(hFig,'materials.png','-dpng','-r300')
    end


    clear A
    clear B
    clear C
    clear D
    clear closeStatus
    clear convFactor
    clear fid
    clear hFig
    clear hVec1
    clear iMaterial
    clear legStrs
    clear m
    clear nMaterials
    clear parAlpha
    clear parN
    clear pos
    clear relCond
    clear seekStatus
    clear thetaAct
    clear thetaRes
    clear thetaSat
    clear r
    clear pos
    clear maxWidth
    
end % if viz



