function PE_ram_run(svp_path,svg_path,out_root,xy_src,xy_end,zs_lims,f_lims,safety_fac,delr,delz,npt,zmplt)
%
% Call PE_ram_run(svp_path,svg_path,out_root,xy_src,xy_end,zs_lims,f_lims,safety_fac,delr,delz,npt,zmplt);
%
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
global SVG BATHY;
%
if nargin<8, safety_fac=1; end
if nargin<9, delr=0; end
if nargin<10, delz=0; end
if nargin<11, npt=0; end
if nargin<12, zmplt=0; end
%
SVP=adi_svp_read(svp_path);
GAP=adi_svp_interpret(SVP);
zs=adi_lim_eval(zs_lims); nzs=length(zs);
freq=adi_lim_eval(f_lims); nfreq=length(freq); freq0=mean(freq);
%
% Read in all SVPs from _svg file used:
svg_open(svg_path);
nprof=fscanf(SVG.fid,'%i',1);
svp_svg=zeros(100,2,nprof); nz_svg=zeros(1,nprof);
for jprof=1:nprof,
	nz_svg(jprof)=fscanf(SVG.fid,'%i',1);
	zc=fscanf(SVG.fid,'%f',2*nz_svg(jprof));
	svp_svg(1:nz_svg(jprof),1:2,jprof)=reshape(zc,2,nz_svg(jprof))';
end
%
iifig=0; iipe=1; iikm=1;		% iipe=1 ==> don't limit hw to hmax done for modes;
if npt<=0, npt=60; end
[ii_a,r_a,h_a,ip_a]=adi_prov_rad(xy_src,xy_end,iikm,npt,iifig,iipe);
%
% nbathy=1; ksec(1)=1; nmid=4;
% r1=0; h1==interp2(BATHY.xax,BATHY.yax,BATHY.hw,xy_src(1),xy_src(2)); 
% hw_vs_rng(1,:)=[r1 h1]; 
% for jsec=2:length(ii_a),
% 	r2=r_a(jsec); h2=h_a(jsec);
% 	for jmid=1:nmid,
% 		fac=jmid/nmid; rmid=r1 + fac*(r2-r1); hmid=h1 + fac*(h2-h1); 
% 		nbathy=nbathy + 1; hw_vs_rng(nbathy,:)=[1000*rmid,hmid];
% 		ksec(nbathy)=jsec-1;
% 	end
% 	r1=r2; h1=h2;		% Missing before???
% end
x_ax=mylinspace(xy_src(1),xy_end(1),npt)'; y_ax=mylinspace(xy_src(2),xy_end(2),npt)';
r_ax=sqrt((x_ax-x_ax(1)).^2 + (y_ax-y_ax(1)).^2); rmax=r_ax(end); rmax_m=1000*rmax; 
% Interpolate bathymetry along radial:
h_ax=interp2(BATHY.xax,BATHY.yax,BATHY.hw,x_ax,y_ax);
hw_vs_rng=[1000*r_ax,h_ax];
nbathy=npt;
% Fill in array containing TL sections bathy points lie in:
dr_ax=r_ax(2)-r_ax(1);
ndr_ax=floor((r_a+.001)/dr_ax) + 1; j1=1;
nsec=length(~isnan(ii_a));
for jsec=1:nsec-1,
	j2=ndr_ax(jsec+1);
	if j2>=j1, ksec(j1:j2)=jsec; end
	j1=j2+1;
end
if nbathy>j1, ksec(j1:nbathy)=nsec-1; end

%
% Find minimum sound speed in first profile and bottom (use as PE reference sound speed):
jline1=ii_a(1); svp_sec1=svp_svg(1:nz_svg(jline1),:,jline1);
cbot=squeeze(GAP.geo(:,1,:));
cmin=min([svp_sec1(:,2); cbot(:)]);
%
% Compute PE sampling parameters:
lambda=1500/freq0;
nlam_bot=10; nlam_false=3;
hbot=sum(GAP.h);
zmax=max(h_a) + hbot + (nlam_bot+nlam_false)*lambda;
if zmplt==0, zmplt=max(h_a) + hbot + round(10*lambda); end
dr=lambda/safety_fac; dz=lambda/(4*safety_fac);
rmax_m=1000*rmax;
if delr>0, ndr=max(1,round(delr/dr)); else, ndr=max(1,ceil(10*lambda/dr)); end
if delz>0, ndz=max(1,round(delz/dz)); else, ndz=max(1,ceil(.25*lambda/dz)); end
%
ram_root=['ramPE_runs']; ram_infile='ramPE_ram.in';
for jf=1:nfreq,
	for jzs=1:nzs,
        fid1=fopen(ram_infile,'w'); 
		if fid1<0, error(['CANNOT OPEN RAM INPUT FILE ',ram_out,'; CHECK NOT IN USE.']); end
		fprintf(fid1,'BATHY = %s; GDEM = %s; SVG = %s\n',SVG.bathy_root,SVG.gdem_root,SVG.out_root);
		fprintf(fid1,'%9.2f\t%9.2f\t%9.2f\t\t\t freq zs zr\n',[freq(jf) zs(jzs) 0]);
		fprintf(fid1,'%9.2f\t%9.4f\t%9i\t\t\t rmax_m dr ndr\n',[rmax_m dr ndr]);
		fprintf(fid1,'%9.2f\t%9.4f\t%9i\t%9.2f\t\t zmax dz ndz zmplt\n',[zmax dz ndz zmplt]);
		fprintf(fid1,'%9.2f\t%9i\t%9i\t%9.2f\t\t c0 npade ns rs\n',[cmin 5 1 0]);
		fprintf(fid1,'%9.2f %9.2f \t\t\t\t rbathy zbathy \n',hw_vs_rng(1,:));
% Write out bathymetry:
		for jb=2:nbathy, fprintf(fid1,'%9.2f %9.2f \n',hw_vs_rng(jb,:)); end
		fprintf(fid1,'-1  -1 \n');
% Loop over bathymetry points:
		for jbathy=1:nbathy,
			jsec=ksec(jbathy);
			jline1=ii_a(jsec); jline2=ii_a(jsec+1);
			svp_sec1=svp_svg(1:nz_svg(jline1),:,jline1); 
			svp_sec2=svp_svg(1:nz_svg(jline2),:,jline2); 
			rmid=hw_vs_rng(jbathy,1); hmid=hw_vs_rng(jbathy,2);
			if h_a(jsec)>h_a(jsec+1), svp0=svp_interp(hmid,svp_sec1); else, svp0=svp_interp(hmid,svp_sec2); end
			ndepth=size(svp0,1);
			if jbathy>1, fprintf(fid1,'%9.2f \t\t\t\t\t rbathy \n',rmid); end
% Print out profile for this range:
			PE_prof_write(fid1,svp0,ndepth,GAP,lambda,nlam_bot,nlam_false,freq(jf));
		end
		fclose(fid1);
% Run RAM PE:
		ramp(ram_infile,ram_root);
%        
		SD=hdfsd_open([ram_root,'.hdf']);
		[TL,status]=hdfsd_read(SD,[],1);
		[PHASE,status]=hdfsd_read(SD,[],2);
        if jf==1 & jzs==1,
            nzr=TL.dims(1); nrng=TL.dims(2);
            hdf_root=out_root; SDw=[]; jset=1; dlab='PE TL - dB';
            zr_ax=TL.scales{1}; rng_ax=TL.scales{2};
			dims=[nzr,nrng,nzs,nfreq]; scales={zr_ax rng_ax zs(:)' freq(:)'};
			labels={'Depth - m', 'Range - km', 'Source Depth - m','Frequency - Hz'};
			[SDw,stat]=hdfsd_write_open(hdf_root,SDw,1,dlab,dims,scales,labels);
			dlab='PE TL Phase - deg';
			[SDw,stat]=hdfsd_write_open(hdf_root,SDw,2,dlab,dims,scales,labels);
        end
% 		disp(['Writing PE TL file ',hdf_root,'.hdf']);
		[stat]=hdfsd_write(SDw,1,{[] [] [jzs,jzs] [jf jf]},TL.data);
% 		disp(['Finished writing PE TL file ',hdf_root,'.hdf']);
		[stat]=hdfsd_write(SDw,2,{[] [] [jzs,jzs] [jf jf]},PHASE.data);
% 		disp(['Finished writing PE TL phase ',hdf_root,'.hdf']);
 		disp(['Finished source depth = ',int2str(round(zs(jzs))),', f = ',int2str(round(freq(jf)))]);
%         TL_dB(:,:,jzs,jf)=TL.data;
%         TL_phase(:,:,jzs,jf)=PHASE.data;
	end
end
%
hdfsd_write_close(SDw);
%
% eval(['!del ',ram_root,'*']);
%     
%
disp(['Finished writing PE TL file ',hdf_root,'.hdf']);
%
return
%
%%%
%
function svp0=svp_interp(h0,svp);
%
ii=find(~isnan(svp(:,1))); nsvp=ii(length(ii)); hsvp=svp(nsvp,1);
if h0<hsvp,
% Cut SVP at h0:
	ii=find(svp(:,1)>h0); n=ii(1); svp0=svp(1:n,:); nn=[n-1:n];
	svp0(n,2)=interp_1d(svp(nn,1),svp(nn,2),h0);
	svp0(n,1)=h0;
elseif h0>hsvp,
% Extrapolate SVP to h0:
	nn=nsvp-1:nsvp;
	c0=interp_1d(svp(nn,1),svp(nn,2),h0);
	svp0=[svp(1:nsvp,:); h0,c0];
else
	svp0=svp(1:nsvp,:);
end
%
%%%
function ci=interp_1d(z,c,zi)
% Straightforward 1-D interpolation.
dz=max(1e-3,z(2)-z(1));
fac=(c(2)-c(1))/dz;
ci=c(1) + (zi-z(1)).*fac;

