function Pmv10

% Pmv10  - reads static shift values from t?.dat files and applies them
%          in pseudosection plots to both, data and response
%        - app & phas curves can be plotted for both modes together
%        - error corrected for cases where only tip is inverted
%        - plots difference (log) between two models
%        - plots sensitivity map, if available

% Pmv9   organizes reading from parameter file (*.par)

% pmv8   heavily reorganized (ws)
%        - handles all three modes (te/tm/hz)
%        - rsp file (output of rund2inv) is analysed directly
%        - modes can have different (sub)sets of frequency/site
%        - log file is read 
%          => used error floors  (only for conjugate=n not in log file),
%          tau, chisq, rms, roughness... read and the latters 
%          can be plotted.

% pmv7   Error floor applied before plotting as option
%        In pseudosections, NaN's can be filled with next highest frequency
%        This is only applied as far as the last frequency.

% pmv6   Fixed bugs in statics (both he last frquencyTE and TM done)
%        rs for residual pseudosections .... with error floor option
%        Separate options for HZ data in plotdat
%        Added subroutine doplot_hz

% pmv5   plots rho vs depth at each site in turn (in "resistivity model")
%        Colour and shading options in model plot hardwired

% pmv4   latest version of Mackie code (march6 2003) has more decimal
%        places in *rsp file.  Changes made in load_nlcgout

% pmv3    computes static if rho doesn't fit

% This version Feb 2002. Plots with new version of Mackies code that
% gives files names from root.

%=======================================================================
% June 2002 : Plots Profiles at a specific period for MT or H_z data
%=======================================================================

% function plot_nlcg_out
%
% Plots nlcg6 responses, data and model
% 
% If a file scales.m exists in working directory, then you can set values
%  for following:
%  scale   = [minY maxY minZ maxZ]; % (kilometers)
%            (minZ and maxZ are optional and not used by this function)
%  rholims = [minlog10rho maxlog10rho]); % (Ohm-m)
%  phalims = [minpha maxpha]); % degrees
% Otherwise these are set automatically and you can alter the values and
%  save the file for future use. Note that missing or bad data in this
%  file may cause the function to return to defaults or quit with an error.
%
% If a file pal.m exists in working directory, it will be executed to
%  load a custom colormap array "cmap". Otherwise a modified jet is used.
%

disp(' ');
disp('Plots NLCG6 output:');
disp(' ');
disp('pseudosection / sounding curves / resistivity model.');
disp(' ');
disp('requires: NLCG_6_9 data input files (.te/.tm/.tip/.hz) and');
disp('          parameter (.par), model (.mdl), response (.rsp) and log (.log) output files.');
disp(' ');

global statics lstatic
global err_floor 
global cmode type kmode

cmode = {'TM';'TE';'HZ';'TMTE'}; 
type = {'rho','pha';'rho','pha';'rel','img'};
version = '6_9';

[parfile,outpath] = uigetfile({'*.par'},'pick inversion parameter file');
parfile=[outpath,'/',parfile];
if exist(parfile,'file')
    disp(['reading parameter file: ',parfile]);
    [seed,kmode,modfile,datfile,tau,err_floor,statshift,indsitmod] = ...
        load_par(parfile);
else
    error([parfile,' does not exist!'])
end
seedp=[outpath,seed,'_',version];
modfile=[outpath,'/',modfile];

% convert error floors to log10 scale
% drho in log10 space: 
% dlog10(rho) = d(ln(rho))/ln(10) = rel.d(rho)/ln(10)
err_floor(1:2,1) = 0.434.*err_floor(1:2,1)./100;
% rund2inv doc tells me (ws) that phase error is in percent of 
% equivalent ln(rho), and thus: 1/(2*100)*180/pi = 0.286 in deg
err_floor(1:2,2) = .286.*err_floor(1:2,2);
% don't change Hz error floor.

% load the response file which will be called dat
if exist([seedp,'.rsp'],'file')
    disp(['reading response file: ',seedp,'.rsp']);
    [dat,per,indsit,sitrms] = load_nlcgout([seedp,'.rsp'],kmode,indsitmod);
else
    error(['File ',seedp,'.rsp', ' does not exist!'])
end
% get freq
freq = 1./per; logfreq = log10(freq);
nfreq = length(freq)
nsites=length(indsit)
   
% load the model file
if exist(modfile,'file')
    [ny,nz,dy,dz,model] = load_mdl(modfile,version);
else
    error(['File: ',modfile,' does not exist'])
end

% read objfun (rms, roughness...) as a function of iteration
nit=0;objfun=zeros(1,6);
if exist([seedp,'.log'],'file')
    [objfun] = load_log([seedp,'.log'],err_floor);
    nit=size(objfun,1)-1;
else
    warning(['File: ',[seedp,'.log'],' does not exist'])
end

% compute y's and z's at block edges
y=zeros(1,ny+1);
for k=1:ny  y(k+1)=y(k)+dy(k); end
z=zeros(1,nz+1);
for k=1:nz  z(k+1)=z(k)+dz(k); end

% compute y's at sites
ysite=y(indsit)+dy(indsit)'/2;

y0=0;
% sites, y: site location wrt 1st site (+y0) in km
sites=ysite-ysite(1)+y0;
sites=sites/1000;
y=sites;

% Initialize statics
lstatic = zeros(2,nsites);
statics = zeros(2,nsites); 

% set defaults
ymin=min(y);     	ymax=max(y);             dy=ymax-ymin;
ymin=ymin-.1*dy;	ymax=ymax+.1*dy;
rhomin=0;        	rhomax=3;                
phamin=0;        	phamax=90;

% load errors matrix
err = zeros(3,nsites,nfreq,2);
for imode=1:3
    if kmode(imode)==1
        name=[outpath,'/',char(datfile(imode))];
        if exist(name,'file')
            [err]=load_err(err,name,imode,indsit,indsitmod,per);
        else
            warning(['File: ',name,' does not exist'])
        end
    end
end
% modes,freqs,sites,[meas_1 meas_2]
err=permute(err,[1 3 2 4]);

% load static shifts
lstatic=zeros(2,nsites);
statics=zeros(2,nsites);
if statshift == 'y'
    [lstatic,statics]=load_stat(kmode,cmode,nsites,indsit,indsitmod,...
                                outpath)
end

dat=reshape(dat(:,:,1:12),[nsites,nfreq,4,3]);
dat=permute(dat,[4 2 1 3]);
% modes,freqs,sites,[obs_1 obs_2 meas_1 meas_2],mode

% high (>1E+5) errors -> NaN
dat(err>1E+5)=NaN;
err(err>1E+5)=NaN;

% modes are 'found' when there's at least one data point ~= NaN
no_nan = ~isnan(dat);
kmode = sum(sum(sum(no_nan,4),3),2)>0 % e.g. [1 0 1] <= TM and HZ modes inverted
mtmode = find(kmode(1:2)>0);          % e.g. [1] ...

% scale error of rho & pha & tzy
% d(ln(rho)) => d(log10(rho)) = d(ln(rho))/ln(10)
err(mtmode,:,:,1) = 0.43.*err(mtmode,:,:,1);
err(mtmode,:,:,2) = (180/pi).*err(mtmode,:,:,2);
% err(3,:,:,:)   = err(3,:,:,:).*0.707;
% I (ws) guess the error in the input file does not refer to the 
% absolute error of the complex number but its real and imag parts

% log10(rho), pha=-pha
dat(1:2,:,:,[1:2:3]) = log10(dat(1:2,:,:,[1:2:3]));
dat(1:2,:,:,[2:2:4]) = -dat(1:2,:,:,[2:2:4]);
        
% extend data matrices for pcolor
err(:,nfreq+1,:,:) = err(:,nfreq,:,:);
err(:,:,1+nsites,:) = err(:,:,nsites,:);
dat(:,nfreq+1,:,:) = dat(:,nfreq,:,:);
dat(:,:,nsites+1,:) = dat(:,:,nsites,:);

choice=menu('Apply error floor to data?','Yes','No');
%=============================================================================
% Apply error floors
%=============================================================================
if choice==1
    for imode=1:3
        if kmode(imode)
            for itype=1:2
                terr = err(imode,:,:,itype);
                lindex = find(terr<=err_floor(imode,itype));
                hindex = find(~isnan(terr));
                terr(lindex)=err_floor(imode,itype);
                err(imode,:,:,itype)=terr;
                nlow = length(lindex); nhigh = length(hindex);
                if nhigh>0 percent=100*nlow/nhigh; else percent=0; end
                disp(['Percent of ',upper(char(type(imode,itype))),...
                  ' data (',char(cmode(imode)),' mode) below error floor ',...
                   num2str(err_floor(imode,itype)),': ',num2str(percent)]);
            end
        end
    end 
end

hmode=[];
for imode=1:3
    if kmode(imode)>0 hmode=[hmode,' ',char(cmode(imode))]; end
end

% loop for ever
while 1
    choice=menu('plot option','plot pseudosection','plot sounding curves',...
                'profiles','resistivity model','sensitivity map',...
                'objective function','quit');
    switch choice
        case {1,2,3}
            lindex=find(kmode>0);
            if length(lindex) > 1
                % TM and TE mode? => TMTE (together)
                cmmode=cmode(lindex);
                if choice==2 & find(lindex==1) & find(lindex==2) 
                    cmmode{length(lindex)+1}='TMTE';
                end
                imode = menu('choose mode',cmmode);
                if imode==length(lindex)+1
                    imode=4;
                else
                    imode = lindex(imode);
                end
            else
                imode = lindex(1);
            end
            if choice==1
                plotdat(logfreq,sites,dat,err,imode,outpath)
            elseif choice==2
                docurves(logfreq,sites,dat,err,imode,outpath)             
                dostatics(outpath)
            elseif choice==3
                doprofiles(logfreq,sites,dat,err,ymin,ymax,imode,outpath)
            end
        case 4
            header=['model: ',seed,' |',hmode,' | tau=',num2str(tau),...
                    ' | rms=',num2str(objfun(nit+1,2)),' | ',num2str(nit),' iters'];
            %header=''
            nlcg_mdl(seedp,header,outpath)
        case 5
            header=['sens: ',seed,' |',hmode,' | tau=',num2str(tau),...
                         ' | rms=',num2str(objfun(nit+1,2)),' | ',num2str(nit),' iters'];
            %header=''
            nlcg_sns(seedp,header,outpath)
        case 6
            plot_objfun(objfun,outpath)
        case 7
            % close all
            return
    end
end

%==========================================================================
function [seed,kmode,modfile,datfile,tau,err_floor,statshift,indsitmod]=...
              load_par(name)
%==========================================================================
% read inversion parameter file
clear datfile err_floor 
indsitmod=zeros(3,1); err_floor=zeros(3,2);

fid=fopen(name,'r');

line=fgetl(fid); seed=sscanf(line,'%s',1);
line=fgetl(fid); kmode(1)=sscanf(line,'%s',1);
line=fgetl(fid); kmode(2)=sscanf(line,'%s',1);
line=fgetl(fid); kmode(3)=sscanf(line,'%s',1);
kmode=kmode=='y'
line=fgetl(fid); modfile=sscanf(line,'%s',1);
line=fgetl(fid); devmod=sscanf(line,'%d',1);
modind=find(kmode>0);
for imode=1:length(modind)
    line=fgetl(fid); datfile{modind(imode)}=sscanf(line,'%s',1);
end
if kmode(3)==1 line=fgetl(fid); conj=sscanf(line,'%s',1); end
line=fgetl(fid); uni_std_lap=sscanf(line,'%d',1);
line=fgetl(fid); grd_lap_reg=sscanf(line,'%d',1);
line=fgetl(fid); tau=sscanf(line,'%d',1);
line=fgetl(fid); sensout=sscanf(line,'%s',1);
if kmode(1)==1 line=fgetl(fid); err_floor(1,1)=sscanf(line,'%f',1)'; 
               line=fgetl(fid); err_floor(1,2)=sscanf(line,'%f',1)';
end
if kmode(2)==1 line=fgetl(fid); err_floor(2,1)=sscanf(line,'%f',1)'; 
               line=fgetl(fid); err_floor(2,2)=sscanf(line,'%f',1)';
end
if kmode(3)==1 line=fgetl(fid); err_floor(3,1)=sscanf(line,'%f',1); 
                                err_floor(3,2)=err_floor(3,1);
end
line=fgetl(fid); parfix=sscanf(line,'%s',1);
if (parfix=='y') 
    line=fgetl(fid); fixfile=sscanf(line,'%s',1);
    line=fgetl(fid); fixdamp=sscanf(line,'%d',1);
end

if (kmode(1)==1 | kmode(2)==1)
    line=fgetl(fid); statshift=sscanf(line,'%s',1)
    if (statshift=='y') 
        line=fgetl(fid); statvar=sscanf(line,'%d',1);
        line=fgetl(fid); statdamp=sscanf(line,'%d',1);
    end
else
    statshift='n'
end
line=fgetl(fid); maxit=sscanf(line,'%d',1)

imode=0; inodep=1E6;
while 1
    line=fgets(fid);
    if line==-1 break; end
    if strfind(line,'END') break; end
    inode=str2num(line);
    if inode<inodep
        imode=imode+1; 
        isite=1;
    else
        isite=isite+1;
    end
    indsitmod(modind(imode),isite)=inode;
    inodep=inode;
end

fclose(fid);
return

%==========================================================================
function [objfun]=load_log(name,err_floor)
%==========================================================================
% get used error_floor and objective function
% objf(iter,1:6)
% 1            2                         3       4                5               6
% s1 [=chisq], rms [=sqrt(chisq/ndpts)], s2/tau, s2 [=roughness], s3[=closeness], s [=objfunc]

fid=fopen(name,'r');
while 1
    line=fgetl(fid);
    if line==-1 break; end
    word=strread(line,'%s');
    if strfind(char(word(1)),'iter')
        while 1
            line=fgetl(fid);
            if line==-1 fclose(fid); return; end
            [temp]=sscanf(line,'%f');
            objfun(temp(1)+1,1:6)=temp(2:7)';
        end
    end
end
        
%==========================================================================
function [ny,nz,dy,dz,rho]=load_mdl(name,version)
%==========================================================================
fid=fopen(name,'r');

A=fscanf(fid,'%i',2);
ny=A(1);
nz=A(2);
dy=fscanf(fid,'%f',ny);
dz=fscanf(fid,'%f',nz);
if version==6
   fscanf(fid,'%i',1);
end
rho=fscanf(fid,'%f',ny*nz);
rho=reshape(rho,ny,nz)';
rho=log10(rho);

fclose(fid);
return

%==========================================================================
function [ny,nz,dy,dz,sens]=load_sns(name,version)
%==========================================================================
fid=fopen(name,'r');

A=fscanf(fid,'%i',2);
ny=A(1);
nz=A(2);
dy=fscanf(fid,'%f',ny);
dz=fscanf(fid,'%f',nz);
if version==6
   fscanf(fid,'%i',1);
end
sens=fscanf(fid,'%f',ny*nz);
sens=reshape(sens,ny,nz)';
sens=log10(sens);

fclose(fid);
return

%=============================================================================
function [dat,per,indsit,sitrms]=load_nlcgout(name,kmode,indsitmod)
%=============================================================================
% this function reads nlcg output.data into a matrix dat;
% dat matrix has 13 colomuns :
%       obsTMrho  obsTMphs    calTMrho  calTMphs ...
%   ... obsTErho  obsTEphs    calTErho  calTEphs ...
%   ... obsHzreal obHzimag    calHzreal calHzimag     Period
% there will be freq*sites rows

fid=fopen(name,'r');
line=1; isite=0;
while line ~= -1
    line=fgetl(fid);
    if findstr(line,'block')
        isite=isite+1;
        word=strread(line,'%s');
        indsit(isite)=str2num(char(word(4)));
        line=fgetl(fid); line=fgetl(fid);
        word=strread(line,'%s');
        sitrms(isite)=str2num(char(word(7)));
    end
    if findstr(line,'===')
        iper=0;
        % which modes for this site (indsitmod from parameter file)
      %smode=sum(indsitmod==indsit(isite),2)>0;
        % results in indices:
      %dind=[reshape((smode*ones(1,4).*[1:4; 5:8; 9:12])',1,[]) 13];
        % only indices>0:
        %dind=find(dind>0);
        % start of data section of site isite
        while 1
            line=fgetl(fid);
            if ~(length(line)>30) break; end
            iper=iper+1;
            word=strread(line,'%s');
            for iword=1:length(word)
                if findstr(char(word(iword)),'*') word{iword}='NaN'; end
            end
      %dat(isite,iper,dind)=str2num(char(word(dind)));
            dat(isite,iper,:)=str2num(char(word));
        end
    end
end

% period
per=sort(reshape(dat(:,:,13),1,[]));
per=per(find(per>0));
per=[per(diff(per)>0), per(length(per))]

% adjust data indices according to period (not sure if ever needed)
for isite=1:size(dat,1)
    tper=squeeze(dat(isite,:,13));
    dat(isite,:,13)=per;
    tdat=squeeze(dat(isite,:,1:12));
    % set data to zero and rewrite
    dat(isite,:,1:12)=NaN;
    for iper=1:length(per)
        for itper=1:length(tper)
            if tper(itper)==per(iper)
                dat(isite,iper,1:12)=tdat(itper,:);
            end
        end
    end
end

fclose(fid);
return

%==========================================================================
function [err]=load_err(err,name,imode,indsit,indsitmod,per)
%==========================================================================
% reads nlcg input data:
%
% data error matrix from input data
% rho errors:   ln(rho)
% phase errors: radians
%
fid=fopen(name,'r');

line=1;
isite=0;
while line ~= -1
    line=fgetl(fid);
    k=findstr(line,'.');
    j=length(k);
    if j>2
        if next_site==1
            next_site=0;
            % isite: index for this site
            isite=isite+1; iper=0;
            % ksite: overall site index
            ksite=find(indsit==indsitmod(imode,isite));
            % was data at this station found in rsp file?
        end
        iper=iper+1;
        % get correct period index
        tfreq=str2num(char(strread(line,'%s',1)));
        for iper=iper:length(per)
            if (tfreq/per(iper)-1 < 0.001) break; end
        end
        tchar=line(findstr(line,')')+1:length(line));
        temp = sscanf(tchar,'%f');
        if (imode==3) temp(2)=temp(1); end
        err(imode,ksite,iper,:) = temp;
    else
        next_site=1;
    end
end

fclose(fid);
return

%==========================================================================
function [lstatic,statics]=load_stat(kmode,cmode,nsites,indsit,indsitmod,...
                                     outpath)
% what data are inverted for static shift? 
% read statics.tx ==> lstatic
% read tx.dat     ==> statics
global statics lstatic
kmode
for imode=1:2
    if kmode(imode)==1
        % which mode & station has been inverted for static shift?
        name=[outpath,'statics.',lower(char(cmode(imode)))];
        if exist(name,'file')
            ltstatic=load(name);
            for isite=1:length(ltstatic)
                ksite=find(indsit==indsitmod(imode,isite));
                lstatic(imode,ksite)=ltstatic(isite);
            end
            disp(['Found: ',name]);
        else
            % no statics.tx file found 
            lstatic(imode,1:nsites)=ones(1,nsites);
            disp(['Not found: ',name]);
            disp(['Set lstatic(',imode,',:) to ones']);
        end
        if sum(lstatic(imode,:))>0
            % what are the actual static shift values?
            name=[outpath,lower(char(cmode(imode))),'.dat'];
            if exist(name,'file')
                disp(['Found: ',name]);
                fid=fopen(name,'r');
                line=1;
                isite=0;
                while line ~= -1
                    line=fgetl(fid);
                    k=findstr(line,'.');
                    if length(k)==1
                        isite=isite+1;
                        temp = sscanf(line,'%f');
                        ksite=find(indsit==indsitmod(imode,isite));
                        % statics found are dat(pred)/dat(obs) => log10
                        statics(imode,ksite)=log10(temp(2));
                        %================================================
                        % note that statics is not used here since it is
                        % applied to both, data and response in 
                        % *.rsp, t?_pred.dat, and t?.dat files.
                        %================================================
                        %statics(imode,ksite)=0;
                    end
                end
                fclose(fid);
            else
                disp(['Not found: ',name]);
                disp(['Set statics(',imode,',:) to zeros']);
                statics(imode,ksite)=zeros(nsites);
            end
        end
    end
end
return

%==========================================================================
function [indsit,sitrms]=load_sites(name)
%==========================================================================
% get the number of sites and index

fid=fopen(name,'r');
line=1; 
m=0;
while line ~= -1
    line=fgetl(fid);
    k=findstr(line,'block');
    if ~isempty(k)
        m=m+1;
        indsit(m)=sscanf(line(k+5:max(size(line))),'%i');
        line=fgetl(fid); line=fgetl(fid);
        word=strread(line,'%s');
        sitrms(m)=str2num(char(word(7)));
    end
end
fclose(fid);
return

%==========================================================================
function plotdat(logfreq,sites,dat,err,imode,outpath)
%==========================================================================
% plot Pseudosection part
% dat:  [modes,freq,sites, obs_1 obs_2 cal_1 cal_2]
% err:  [modes,freq,sites, obs_1 obs_2];

global statics lstatic
global err_floor;
global cmode type

% set the interpolation type
choice = menu('Fill NaN values?','Yes','No');
if choice == 1; ifill=1; else; ifill=0; end

nfreq=size(dat,2); 
nsites=size(dat,3);

% copy relevant data in new array
tdat=squeeze(dat(imode,:,:,:));
terr=squeeze(err(imode,:,:,:));

%=============================================================================
% Apply static shifts to measured MT impedance data
%=============================================================================
if imode<3
    for is=1:nsites-1
        tdat(:,is,1) = tdat(:,is,1)-statics(imode,is);
        % static shift free parameter in inversion?
        % yes ==> also shift responses
        if sum(lstatic(imode,:))>0
            tdat(:,is,3) = tdat(:,is,3)-statics(imode,is);
        end
    end
end 

%=============================================================================
% Set shading type
%=============================================================================
while 1
    choice=menu('shading','interpolated','flat','go back to main');
    if choice==1; shade=1; elseif choice==2; shade=2; else; return
end
% modify y and f so that the pcolor boxes 'span' nodes
y=sites;         nsites=max(size(sites));
freq=logfreq;    nfreq=max(size(freq)); 
dy=diff(y);
yy(1)=y(1)-dy(1)/2;    yy(2:nsites)=y(1:nsites-1)+dy/2; 
yy(nsites+1)=y(nsites)+dy(nsites-1)/2;
y=yy;   clear yy, dy;

nsites=nsites+1;
f=freq;  df=diff(f);
ff(1)=f(1)-df(1)/2;
ff(2:nfreq)=f(1:nfreq-1)+df/2;  ff(nfreq+1)=f(nfreq)+df(nfreq-1)/2;
f=ff;      nfreq=nfreq+1;  clear ff, df;

% set limits 
rhomin=0; rhomax=3;     phamin=0; phamax=90;     hzmin=-0.3; hzmax=0.3;
colbw=0;

% Control axis in pcolor
plot19 = [min(y),max(y),min(f),max(f)];
resid_plot = [-5,5];
done=0;

while ~done
    lims=[-rhomax -rhomin; phamin phamax; hzmin hzmax];
    %======================================================================
    % Plot data and responses
    %======================================================================
    h1 = figure;
    
    if colbw==0
        % define colormap (old!)
        % colmap=menu('color map','John Color','Matlab Jet Color')
        colmap=2;
        if colmap==1;              
            color=[  0    0    0    0  154  255  255  255  255;...
                     0  191  255  255  205  255  220  165    0;...
                   255  255  255    0   50    0    0    0    0]';
            cmap=color/255;
        else
            cmap=jet(64);
            %   cmap(20,:)=[];  cmap(19,:)=[];  cmap(13,:)=[]; 
            %   cmap(7,:)=[];   cmap(1,:)=[];   cmap(1,:)=[];   
        end
    else
        cmap=gray(13);
    end
    colormap(cmap);
        
    % loop over types, e.g.: [obs_rho, obs_pha, cal_rho, cal_pha] 
    type_match = [1 3 2 4];
    for itype=1:4
%        if rem(itype,2)==0 

        ktype=type_match(itype);
        subplot(4,1,itype)
        %subplot(2,1,itype/2) 

        % Fill NaN values
        if ifill==1 & ktype<3
            % Determine the last frequency w/ data
            for is=1:nsites
                if sum(tdat(:,is,ktype))>0
                    no_data=1;
                    nlast(is)=max(find(tdat(:,is,ktype)>0));
                else
                    no_data=0;
                    nlast(is)=nfreq;
                end
                for ifreq=1:nlast(is)
                    if terr(ifreq,is,ktype) >= 1E+5;
                        if ifreq>1 & no_data==0
                            tdat(ifreq,is,ktype) = tdat(ifreq-1,is,ktype);  
                        elseif is>1 | no_data==1
                            tdat(ifreq,is,ktype) = tdat(ifreq,is-1,ktype);  
                        end
                    end
                end    
            end      
        end
    
        if (itype<3 & imode<3) ksign=-1; else ksign=1; end
        
        pcolor(y, f, ksign.*tdat(:,:,ktype));
        axis(plot19);
        hold on
    
        plot(sites,f(1)*ones(nsites-1,1),'kv');
        if imode<3 & itype<3   klims=1;  elseif ...
           imode<3 & itype>2   klims=2;  elseif ....
           imode==3            klims=3;  
        end
        caxis(lims(klims,:)); 
        
        if itype==1 title([char(cmode(imode)),' Mode']); end
    
        colorbar('vert');
        if shade==1; shading interp; else; shading flat; end
        hold off
        % end
    end

    %====================================================
    % Plot residuals
    %====================================================   
    h2 = figure;

    % here: type=[rho abs] or [rel img].
    for itype=1:2
        % subplot(211)
        subplot(2,1,itype)
        temp = (tdat(:,:,itype)-tdat(:,:,itype+2))./terr(:,:,itype);
        if colbw==1 temp=abs(temp); end
        
        pcolor(y, f, temp);
    
        if colbw==0 colormap(cmap); else colormap(flipud(cmap)); end
        
        axis(plot19);
        hold on
        plot(sites,f(2)*ones(nsites-1,1),'kv');
        title([char(cmode(imode)),' mode residuals: ',...
               upper(char(type(imode,itype)))]);
       
        if colbw==0 caxis(resid_plot); else caxis([0 resid_plot(2)]); end;

        colorbar('vert')
        if shade==1; shading interp; else; shading flat; end
        hold off
    end
 
    % change scales
    if imode<3
        choice=menu('Change:','resistivity limits','phase limits','color/bw','nothing');
        if choice==1
            set_params1=1;
            rhomin=input('Enter min log10(resistivity (Ohm-m)): ');
            while 1
                rhomax=input('Enter max log10(resistivity) > min: ');
                if rhomax>rhomin  break; end
            end
        elseif choice==2
            set_params=1;
            phamin=input('Enter min phase (deg): ');
            while 1
                phamax=input('Enter max phase > min phase: ');
                if phamax>phamin  break; end
            end
        elseif choice==3
            colbw=rem(colbw+1,2);
        else
            done=1;
        end
    else
        choice=menu('Change:','tipper limits','color/bw','nothing');
        if choice==1
            set_params1=1;
            hzmin=input('Enter min tipper: ');
            while 1
                hzmax=input('Enter max tipper > min: ');
                if hzmax>hzmin  break; end
            end
        elseif choice==2
            % change: color<->black and white
            colbw=rem(colbw+1,2);
        else
            done=1;
        end
    end    
end

scolbw={'color','bw'};
% print to file
choice=menu('Make Postscript File?','yes','no');
if choice==1
       disp(['Making PS files for ',char(cmode(imode)),' response']);
       print(h1,'-dpsc',[outpath,'/','nlcg_rsp_',lower(char(cmode(imode))),'_',char(scolbw(colbw+1)),'.ps']);  
       print(h2,'-dpsc',[outpath,'/','nlcg_rsp_',lower(char(cmode(imode))),'_',char(scolbw(colbw+1)),'_err.ps']);
end

% close all
return
end

%==========================================================================
function docurves(lfreq,sites,dat,err,imode,outpath)
% plot sounding curves:     
% -> calls: doplot
% dat:  [modes,freq,sites, obs_1 obs_2 cal_1 cal_2]
% err:  [modes,freq,sites, obs_1 obs_2];
global cmode

% loop for ever;
while 1
    choice=menu('plot option','plot all','plot single site','go back to main');
    if choice==1
        plotall=1;
    elseif choice==2
        plotall=0;
    else
        return
    end
    nsites=length(sites);
    % figure out how many sites need to be plotted
    if plotall==1
        kstart=1;
        kend=nsites;
    else
        while 1
            isite=input('Enter site number to plot (0 to exit): ');
            if ~(isite>nsites) break; end
             fprintf(2,'%s\n','Invalid entry: isite > nsites')
        end
        if isite==0 return; end;
        kstart=isite;
        kend=isite;
    end

    nfreq=length(lfreq);
    % plot curves
    for isite=kstart:kend
        figure(20)
        if (imode<4)
            doplot(isite,lfreq,sites,dat,err,imode,outpath);
        else % tm & te mode in one plot (~imode==4)
            doplot(isite,lfreq,sites,dat,err,1,outpath);
            figure(20)
            doplot(isite,lfreq,sites,dat,err,2,outpath);
        end
        eval(['print -dps ',outpath,'/','site_',num2str(isite),'_',...
             lower(char(cmode(imode)))]);
        pause
        close(20);
    end 
end
return


%=============================================================================
function doplot(isite,lfreq,sites,dat,err,imode,outpath)
%=============================================================================
% dat:  [modes,freq,sites, obs_1 obs_2 cal_1 cal_2]
% err:  [modes,freq,sites, obs_1 obs_2];
global statics lstatic
global cmode

nfreq=length(lfreq);

% note that dat is log10(rho)
tdat=squeeze(dat(imode,1:nfreq,isite,:));
terr=squeeze(err(imode,1:nfreq,isite,:));

fmax=ceil(max(lfreq));
fmin=floor(min(lfreq));

sym={'b.';'r.';'g.'}; symA={'b-';'r-';'g-'};
    
h1=subplot(211);

% Compute static shift for TE and TM mode
if imode<3 
    if sum(lstatic(imode,:))==0 
        % Compute static shift
        % exp(mean(ln(observed)-ln(predicted)))
        stat = tdat(:,1)-tdat(:,3);
        lindex = find(stat>=-5 & stat<=5);
        %statics(imode,isite)=mean(stat(lindex));
        statics(imode,isite)=0;
    end 
    % plotted is shifted data
    errorbar(lfreq,tdat(:,1)-statics(imode,isite),...
             terr(:,1),terr(:,1),char(sym(imode)));    
else
    errorbar(lfreq,tdat(:,1),terr(:,1),terr(:,1),char(sym(imode)));
end
hold on
plot(lfreq,tdat(:,3),char(symA(imode)));
xlabel('Log10 Frequency');
if imode~=3
    ylabel('Log10 Resistivity');
    axis([fmin fmax -1 4]);
else
    ylabel('Tzxr');
    axis([fmin fmax -0.75 0.75]);
end
title([' Site # ',num2str(isite),' at: ',num2str(sites(isite)),' km ',...
        char(cmode(imode)),' Mode']);
set(gca,'XDir','reverse','YGrid','on');
%hold off

h2=subplot(212);
    
errorbar(lfreq,tdat(:,2),terr(:,2),terr(:,2),char(sym(imode)));
hold on
plot(lfreq,tdat(:,4),char(symA(imode)));
xlabel('Log10 Frequency');
if imode<3
    ylabel('Phase');
    axis([fmin fmax  0 90]);
else
    ylabel('Tzxi');
    axis([fmin fmax -0.75 0.75]);
end
set(gca,'XDir','reverse','YGrid','on');
%hold off
eval(['print -dps ',outpath,'/','site_',num2str(isite),'_',lower(char(cmode(imode)))]);

return 
 
%=======================================================================================================
function doprofiles(lfreq,sites,dat,err,ymin,ymax,imode,outpath)
% plot profiles
% dat:  [modes,freq,sites, obs_1 obs_2 cal_1 cal_2]
% err:  [modes,freq,sites, obs_1 obs_2];
% loop for ever;
global cmode

nsites=max(size(sites));
iper=input('Enter period to plot (0 to exit): ');
nfreq=max(size(lfreq));

period = 10^(-lfreq(iper))
tdat = squeeze(dat(imode,iper,1:nsites,:));
terr = squeeze(err(imode,iper,1:nsites,:));
 
if imode<3
    tdat(:,1:2:3) = 10*tdat(:,1:2:3);
    terr(:,1) = 2*tdat(:,1).*terr(:,1);
end
 
subplot(2,1,1)

semilogy(sites,tdat(:,3));
hold on
errorbar(sites,tdat(:,1),terr(:,1),'o')
axis([ymin ymax 0.3 1000])
xlabel('distance(km)')
if imode<3
    ylabel([char(cmode(imode)),' apparent resistivity']);
else
    ylabel([char(cmode(imode)),' real part']);
end
    
subplot(2,1,2)

errorbar(sites,tdat(:,2),terr(:,2),'o')
axis([ymin ymax 0 90])
hold on
plot(sites,tdat(:,4))
xlabel('distance(km)')
if imode<3
    ylabel([char(cmode(imode)),' phase (degrees)']);
else
    ylabel([char(cmode(imode)),' imaginary part']);
end

%set(gca,'YGrid','on');
eval(['print -dpsc ',outpath,'/','profile_f',num2str(iper),'_',lower(char(cmode(imode))),'.ps']);
pause

%close all;
return

%=======================================================================================================
function dostatics(outpath)

%Plot static shift as a function of station number
global statics
global kmode

choice=menu('plot option','plot statics','go back to main');
if choice==1
    nsites = size(statics,2);
    ttitle=[];
    if kmode(1)==1   
        plot([1:nsites],statics(1,:),'+');   
        ttitle=[ttitle,'  + = TM mode  '];
        hold on
    end
    if kmode(2)==1   
        plot([1:nsites],statics(2,:),'o'); 
        ttitle=[ttitle,'  o = TE mode  '];
        hold on
    end
    
    figure
    xzero = [0,nsites+1];    yzero = [0,0];
    plot (xzero,yzero,'k:');
    
    ylabel('static');
    xlabel('site number');
    title(ttitle);
    axis([0 nsites -2. 2.0]);
  
    eval(['print -dpsc ',outpath,'/','statics.ps']);
  
    pause
    % close all;  
end
return

%=============================================================================
function nlcg_mdl(outroot,header,outpath)
%=============================================================================
%
%  WARNING: TOPOGRAPHY NOT IMPLEMENTED!!!!
% 
% INPUTS:
%  header is a title for the model plot
%  looks for files: *.mdl and *.rsp
%
%  If a file scales.m exists in working directory, then you can set values
%  for following:
%
%  scale   = [minY maxY minZ maxZ] (kilometers);
%  rholims = [minlog10rho maxlog10rho]) (Ohm-m)
%  phalims
%  VE (vertical exaggeration)
%
%  Otherwise these are set automatically and you can alter the values and
%  save the file for future use. Note that missing or bad data in this
%  file may cause the function to return to defaults or quit with an error.
%  
%  If a file pal.m exists in working directory, it will be executed to
%  load a custom colormap "cmp". Otherwise a modified jet is used.
%

choice=menu('vertical scale','linear','log','go back main');
if choice==1
    zaxis='linear';
elseif choice==2
    zaxis='log';
else 
    % close all
    return
end

version=6;
%header='NLCG Model';

if exist('scales.m','file')
  scales;
  ymin=scale(1);       ymax=scale(2);
  zmin=scale(3);       zmax=scale(4);
  rhomin=rholims(1);   rhomax=rholims(2);
  phamin=phalims(1);   phamax=phalims(2);
  % VE is also set in file scales.m!!
  set_params=0;
else
  set_params=1;
end

% load the model file
name=[outroot,'.mdl']
    
if exist(name,'file')
    [ny,nz,dy,dz,model]=load_mdl(name,version);
else
    error(['File: ',name,' does not exist'])
end

% load the another model file?
lmod2=menu('compare w/ second model?','yes','no');
if lmod2==1
    choice=menu('same first model?','same','other');
    if choice==2
        [mod1fil,mod1path] = uigetfile('*.mdl;*.mod','pick first model');
        mod1file=[mod1path,'/',mod1fil];
        disp(['reading first model file: ',mod1file]);
        [ny,nz,dy,dz,model]=load_mdl(mod1file,version);
    end
    [mod2fil,mod2path] = uigetfile('*.mdl;*.mod','pick second model');
    mod2file=[mod2path,'/',mod2fil];
    disp(['reading second model file: ',mod2file]);
    [ny2,nz2,dy2,dz2,model2]=load_mdl(mod2file,version);
    % relative change
    model=model-model2;
    header=['log[m_1] - log[m_2]'];
    rhomin=-1; rhomax=+1;
end
    
dy=dy/1000; dz=dz/1000; y0=0;  y0=y0/1000;

% load the site indices
name=[outroot,'.rsp'];
if exist(name,'file')
  [indsit,sitrms]=load_sites(name);
else
  error(['File ',name,' does not exist.'])
end

nsites=max(size(indsit));

% compute y's and z's at block edges
y=zeros(1,ny+1);
for k=1:ny  y(k+1)=y(k)+dy(k);   end
z=zeros(1,nz+1);
for k=1:nz  z(k+1)=z(k)+dz(k);   end

% compute y's at sites
ysite=y(indsit)+dy(indsit)'/2;

% offset site(1) to y0
y=y-ysite(1)+y0;
ysite=ysite-ysite(1)+y0;

% defaults
if set_params
     ymax=max(ysite);        ymin=min(ysite); dely=ymax-ymin;
     ymin=ymin-1*dely;       ymax=ymax+1*dely;
%     zmin=min(z);            zmax=max(z);
     zmin=min(z);            zmax=10;
     rhomin=0;               rhomax=3;
     phamin=0;               phamax=90;   phalims=[phamin,phamax];
     % vertical exageration always starts at 1:1
     VE=1;
     % different limits for display of relative differences
     if lmod2==1 
         rhomin=-1; rhomax=+1;
     end
end

% adjust colorbar toggle
adjclb = 0;
colbw = 0;
%choice=menu('plot type','flat','interpolated');
% pmv5 : hardwired
choice=1;
if choice==1  typ='flat';   else  typ='interp';   end

% extend model so that pcolor will plot all elements
model(nz+1,:)=model(nz,:);   model(:,ny+1)=model(:,ny);

% plot the model
done=0;
h1 = figure;
while ~done
    figure(h1);
    clf
    scale=[ymin ymax zmin zmax];
    rholims=[rhomin rhomax];
  
    % specify color map
    % colmap=menu('color map','John Color','Matlab Jet Color');
    colmap=2;
    if colbw==0
        if colmap==1;              
            color=[  0    0    0    0  154  255  255  255  255;...
                     0  191  255  255  205  255  220  165    0;...
                   255  255  255    0   50    0    0    0    0]';
            cmap=color/255;
        else
            cmap=jet(20);
            cmap(20,:) =[];  cmap(19,:)=[];  cmap(13,:)=[]; 
            cmap(7,:)  =[];  cmap(1,:) =[];
            % below: Volcan's favourite
            % cmap=jet(24);
            % cmap(23,:)=[];    cmap(21,:)=[];    cmap(19,:)=[];    cmap(17,:)=[];
            % cmap(15,:)=[];    cmap(13,:)=[];    cmap(11,:)=[];    cmap(9,:)=[];
            %cmap(7,:)=[];     cmap(5,:)=[];     cmap(3,:)=[];
            cmap=jet(24);
            cmap(23,:)=[];    cmap(21,:)=[];    cmap(19,:)=[];    cmap(17,:)=[];
            cmap(15,:)=[];    cmap(13,:)=[];    cmap(11,:)=[];    cmap(9,:)=[];
            cmap(7,:)=[];     cmap(5,:)=[];     cmap(3,:)=[];
        end
    else
        cmap=gray(13);
    end 

    pcolor(y,z,model);
  
    hold on
    caxis(rholims);
    colormap(flipud(cmap));
    if strcmp(zaxis,'log')
        plot(ysite,1.0000001*z(2)*ones(nsites,1),'kv')
    else
        plot(ysite,z(1)*ones(nsites,1),'kv')
    end
    title(header,'Interpreter','none');
    axis(scale);
    xlabel('km');ylabel('km');
    %set(gca,'YTick',[0:50:350])
    hold off
  
    if strcmp(zaxis,'log'); axis square; else; set(gca,'DataAspectRatio',[VE 1 1]); end

    if strcmp(typ,'flat'); shading flat; else; shading interp;  end

    if strcmp(zaxis,'log')
        set(gca,'YScale','log','YDir','reverse')
    else
        set(gca,'YScale','linear','YDir','reverse')
    end

    hcb = colorbar('ver');
    pltpos = get(gca,'Position');
  
    if adjclb==1
        pltpos = get(gca,'Position');
        pratio = get(gca,'PlotBoxAspectRatio');
        pratio = pratio(2)/pratio(1)*1.25;
        clbpos = get(hcb,'Position');
        clbpos(2) = clbpos(2)+(1-pratio)/2*clbpos(4);
        clbpos(4) = pratio*clbpos(4);
        if pratio<1
            set(hcb,'Position',[clbpos(1),clbpos(2),clbpos(3),clbpos(4)]);
            % needed for Matlab7: reset the model plot to old position
            set(gca,'Position',pltpos);            
        end
    end
    if strcmp(zaxis,'linear')
        choice=menu('Change:','resistivity limits','y limits','z limits',...
                    'vertical exageration','adjust colorbar','color/bw','nothing');
    else
        choice=menu('Change:','resistivity limits','y limits',...
                    'z limits','nothing');
    end
    if choice==1
        set_params=1;
        rhomin=input('Enter min log10(resistivity (Ohm-m)): ');
        while 1
            rhomax=input('Enter max log10(resistivity) > min: ');
            if rhomax>rhomin break; end
        end
    elseif choice==2
        set_params=1;
        ymin=input('Enter left edge (km): ');
        while 1
            ymax=input('Enter right edge (> left edge) (km): ');
            if ymax>ymin break; end
        end
    elseif choice==3
        set_params=1;
        while 1
            zmin= input('Enter top edge (>=0) (km): ');
            if zmin>=0 break; end
        end
        while 1
            zmax= input('Enter bottom edge (> top edge) (km): ');
            if zmax>zmin break; end
        end
    elseif choice==4
        if strcmp(zaxis,'linear')
            kk=menu('Vertical Exageration:','1:1','x2','x3','x10');
            if kk==1;       VE=1;
            elseif kk==2;   VE=2*VE;
            elseif kk==3;   VE=3*VE;
            elseif kk==4;   VE=10*VE;
            end    
        else
            done=1;
        end 
    elseif choice==5
        % adjust colorbar toggle
        adjclb=rem(adjclb+1,2);
    elseif choice==6
        % change: color<->black and white
        colbw=rem(colbw+1,2);
    else
        done=1;
    end
end

scolbw={'color','bw'};
%choice=menu('Make Postscript File?','yes','no');
if strcmp(zaxis,'log') 
    modstr='model_LZ_'; 
else
    modstr='model_Z_';
end
choice =1;
if lmod2==1 modstr=['d',modstr]; end
if choice==1
    eval(['print -dpsc ',outpath,'/',modstr,char(scolbw(colbw+1)),'.ps']);
end

zp=menu('Make rho(z)plot?','yes','no');
if zp==1
    choice=menu('All sites?','yes','no');
    while 1
        if choice == 1;
            ismin=1; ismax=nsites;
        else
            ns=input('Enter station number  (0 to end): ');
            ismin=ns; ismax=ns;
            if ns==0 break; end
        end
        for is=ismin:ismax
            h2 = figure;
            rhoz = model(:,indsit(is))';
            plot(z,rhoz,'o-'); axis([0,100,0,3]); title([' Site: ',num2str(is)]);
            xlabel('depth [km]');
            ylabel('log10(rho)');
            
            % Output   z and log10 (rho(z)) to file
            str=['rhoz_',num2str(is),'.dat'];
            fid1=fopen(str,'w+');
            fprintf(fid1,'%9.3f %9.3f \n',[z;rhoz]);
            fclose(fid1);
            
            str=[outpath,'/','rhoz_',num2str(is),'.ps'];
            if choice==2 print(h2,'-dpsc',str); end
            
            pause
        end
        if choice==1 break; end
    end
end

% close all
return

%=============================================================================
function nlcg_sns(outroot,header,outpath)
%=============================================================================
%
% -> basically a copy of NLCG_MDL!
%
%  WARNING: TOPOGRAPHY NOT IMPLEMENTED!!!!
% 
% INPUTS:
%  header is a title for the model plot
%  looks for files: *.mdl, *.rsp, and *.sns
%
%  If a file scales.m exists in working directory, then you can set values
%  for following:
%
%  scale   = [minY maxY minZ maxZ] (kilometers);
%  rholims = [minlog10rho maxlog10rho]) (Ohm-m)
%  phalims
%  VE (vertical exaggeration)
%
%  Otherwise these are set automatically and you can alter the values and
%  save the file for future use. Note that missing or bad data in this
%  file may cause the function to return to defaults or quit with an error.
%  
%  If a file pal.m exists in working directory, it will be executed to
%  load a custom colormap "cmp". Otherwise a modified jet is used.
%

choice=menu('vertical scale','linear','log','go back main');
if choice==1
    zaxis='linear';
elseif choice==2
    zaxis='log';
else 
    % close all
    return
end

version=6;
%header='NLCG Sensitivities';

if exist('scales.m','file')
  scales;
  ymin=scale(1);       ymax=scale(2);
  zmin=scale(3);       zmax=scale(4);
  rhomin=rholims(1);   rhomax=rholims(2);
  phamin=phalims(1);   phamax=phalims(2);
  % VE is also set in file scales.m!!
  set_params=0;
else
  set_params=1;
end

% load the sensitivity file
name=[outroot,'.sns']
if ~exist(name,'file')
    error(['File: ',name,' does not exist'])
end
[ny,nz,dy,dz,sens]=load_sns(name,version);
dy=dy/1000; dz=dz/1000; y0=0;  y0=y0/1000;

% load the site indices
name=[outroot,'.rsp'];
if exist(name,'file')
    [indsit,sitrms]=load_sites(name);
else
    error(['File ',name,' does not exist.'])
end

nsites=max(size(indsit));

% compute y's and z's at block edges
y=zeros(1,ny+1);
for k=1:ny  y(k+1)=y(k)+dy(k);   end
z=zeros(1,nz+1);
for k=1:nz  z(k+1)=z(k)+dz(k);   end

% compute y's at sites
ysite=y(indsit)+dy(indsit)'/2;

% offset site(1) to y0
y=y-ysite(1)+y0;
ysite=ysite-ysite(1)+y0;

% defaults
if set_params
     ymax=max(ysite);        ymin=min(ysite); dely=ymax-ymin;
     ymin=ymin-1*dely;       ymax=ymax+1*dely;
%     zmin=min(z);            zmax=max(z);
     zmin=min(z);            zmax=10;
     % vertical exageration always starts at 1:1
     VE=1;
end
sensmin=-4; sensmax=1;

% adjust colorbar toggle
adjclb = 0;
colbw = 0;
%choice=menu('plot type','flat','interpolated');
% pmv5 : hardwired
choice=1;
if choice==1  typ='flat';   else  typ='interp';   end

% extend model so that pcolor will plot all elements
sens(nz+1,:)=sens(nz,:);   sens(:,ny+1)=sens(:,ny);

% plot the model
done=0;
h1 = figure;
while ~done
    figure(h1);
    clf
    scale=[ymin ymax zmin zmax];
    senslims=[sensmin sensmax];
  
    % specify color map
    % colmap=menu('color map','John Color','Matlab Jet Color');
    colmap=2;
    if colbw==0
        if colmap==1;              
            color=[  0    0    0    0  154  255  255  255  255;...
                     0  191  255  255  205  255  220  165    0;...
                   255  255  255    0   50    0    0    0    0]';
            cmap=color/255;
        else
            cmap=jet(20);
            cmap(20,:) =[];  cmap(19,:)=[];  cmap(13,:)=[]; 
            cmap(7,:)  =[];  cmap(1,:) =[];
            % below: Volcan's favourite
            % cmap=jet(24);
            % cmap(23,:)=[];    cmap(21,:)=[];    cmap(19,:)=[];    cmap(17,:)=[];
            % cmap(15,:)=[];    cmap(13,:)=[];    cmap(11,:)=[];    cmap(9,:)=[];
            %cmap(7,:)=[];     cmap(5,:)=[];     cmap(3,:)=[];
            cmap=jet(24);
            cmap(23,:)=[];    cmap(21,:)=[];    cmap(19,:)=[];    cmap(17,:)=[];
            cmap(15,:)=[];    cmap(13,:)=[];    cmap(11,:)=[];    cmap(9,:)=[];
            cmap(7,:)=[];     cmap(5,:)=[];     cmap(3,:)=[];
        end
    else
        cmap=gray(13);
    end 

    pcolor(y,z,sens);
  
    hold on
    caxis(senslims);
    colormap(cmap);
    if strcmp(zaxis,'log')
        plot(ysite,1.0000001*z(2)*ones(nsites,1),'kv')
    else
        plot(ysite,z(1)*ones(nsites,1),'kv')
    end
    title(header,'Interpreter','none');
    axis(scale);
    xlabel('km');ylabel('km');
    %set(gca,'YTick',[0:50:350])
    hold off
  
    if strcmp(zaxis,'log'); axis square; else; set(gca,'DataAspectRatio',[VE 1 1]); end

    if strcmp(typ,'flat'); shading flat; else; shading interp;  end

    if strcmp(zaxis,'log')
        set(gca,'YScale','log','YDir','reverse')
    else
        set(gca,'YScale','linear','YDir','reverse')
    end

    hcb = colorbar('ver');
    pltpos = get(gca,'Position');
  
    if adjclb==1
        pltpos = get(gca,'Position');
        pratio = get(gca,'PlotBoxAspectRatio');
        pratio = pratio(2)/pratio(1)*1.25;
        clbpos = get(hcb,'Position');
        clbpos(2) = clbpos(2)+(1-pratio)/2*clbpos(4);
        clbpos(4) = pratio*clbpos(4);
        if pratio<1
            set(hcb,'Position',[clbpos(1),clbpos(2),clbpos(3),clbpos(4)]);
            % needed for Matlab7: reset the model plot to old position
            set(gca,'Position',pltpos);            
        end
    end
    if strcmp(zaxis,'linear')
        choice=menu('Change:','sensitivity limits','y limits','z limits',...
                    'vertical exageration','adjust colorbar','color/bw','nothing');
    else
        choice=menu('Change:','sensitivity limits','y limits',...
                    'z limits','nothing');
    end
    if choice==1
        set_params=1;
        sensmin=input('Enter min log10(sensitivity): ');
        while 1
            sensmax=input('Enter max log10(sensitivity) > min: ');
            if sensmax>sensmin break; end
        end
    elseif choice==2
        set_params=1;
        ymin=input('Enter left edge (km): ');
        while 1
            ymax=input('Enter right edge (> left edge) (km): ');
            if ymax>ymin break; end
        end
    elseif choice==3
        set_params=1;
        while 1
            zmin= input('Enter top edge (>=0) (km): ');
            if zmin>=0 break; end
        end
        while 1
            zmax= input('Enter bottom edge (> top edge) (km): ');
            if zmax>zmin break; end
        end
    elseif choice==4
        if strcmp(zaxis,'linear')
            kk=menu('Vertical Exageration:','1:1','x2','x3','x10');
            if kk==1;       VE=1;
            elseif kk==2;   VE=2*VE;
            elseif kk==3;   VE=3*VE;
            elseif kk==4;   VE=10*VE;
            end    
        else
            done=1;
        end 
    elseif choice==5
        % adjust colorbar toggle
        adjclb=rem(adjclb+1,2);
    elseif choice==6
        % change: color<->black and white
        colbw=rem(colbw+1,2);
    else
        done=1;
    end
end

scolbw={'color','bw'};
%choice=menu('Make Postscript File?','yes','no');
if strcmp(zaxis,'log') 
    sensstr='sens_LZ_'; 
else
    sensstr='sens_Z_';
end
choice =1;
if choice==1
    eval(['print -dpsc ',outpath,'/',sensstr,char(scolbw(colbw+1)),'.ps']);
end

zp=menu('Make sens(z)plot?','yes','no');
if zp==1
    choice=menu('All sites?','yes','no');
    while 1
        if choice == 1;
            ismin=1; ismax=nsites;
        else
            ns=input('Enter station number  (0 to end): ');
            ismin=ns; ismax=ns;
            if ns==0 break; end
        end
        for is=ismin:ismax
            h2 = figure;
            sensz = sens(:,indsit(is))';
            plot(z,sensz,'o-'); axis([0,100,0,3]); title([' Site: ',num2str(is)]);
            xlabel('depth [km]');
            ylabel('log10(sensitivity)');
            
            % Output   z and log10 (rho(z)) to file
            str=['sensz_',num2str(is),'.dat'];
            fid1=fopen(str,'w+');
            fprintf(fid1,'%9.3f %9.3f \n',[z;sensz]);
            fclose(fid1);
            
            str=[outpath,'/','sensz_',num2str(is),'.ps'];
            if choice==2 print(h2,'-dpsc',str); end
            
            pause
        end
        if choice==1 break; end
    end
end

% close all
return

%==========================================================================
function plot_objfun(objfun,outpath)
%==========================================================================

nit=size(objfun,1);
cobjfun={'chisq';'rms';'roughness_vs_tau';'roughness';'closeness';...
         'objective_function';'back to main'};
while 1
    choice=menu('Plot',cobjfun);
    if choice==7 return; end
    figure
    plot([0:nit-1],objfun(:,choice));
    title(char(cobjfun(choice)));
    eval(['print -dpsc ',outpath,'/',char(cobjfun(choice)),'.ps']);
end