function PE_convert_ri(ram_root,SVP,GAP,zs,rmax,freq,safety_fac,iires,freq0,delr,delz)
%TL,hw_vs_rng_adi,GAP,jbrg,zs,iiparm,ram_root, ...
%	svp_svg,nz_svg,rmax,zrmax,delz0)
% Call  [hw_vs_rng,svp_src]=PE_convert_ri(TL,hw_vs_rng_adi,GAP,jbrg,zs,iiparm,ram_root, ...
%	svp_svg,nz_svg,rmax,zrmax,delz0);
%
% Produces PE input file for bathymetry and profiles along bearing jbrg
% of the TL structure produced by adi_tl_radial.
%
% freq0 is used to compute lambda so that multiple frequency runs get same output grid;
%
if nargin<9, freq0=freq; end
if nargin<10, delr=0; end
if nargin<11, delz=0; end
%
fid1=fopen([ram_root,'.in'],'w'); 
if fid1<0, 
	error(['USER ERROR: CANNOT OPEN OUTPUT FILE ',ram_root,'; CHECK NOT IN USE.']);
end
%
hw=SVP.ssp(end,1);
% Compute PE sampling parameters:
lambda=1500/freq0;
nlam_bot=10; nlam_false=3;
hbot=sum(GAP.h);
zmax=hw + hbot + (nlam_bot+nlam_false)*lambda; 
zmplt=hw + round(10*lambda);
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
if zs==-999, zs=hw-.01; end
%
% Find minimum sound speed in SSP and bottom:
cbot=squeeze(GAP.geo(:,1,:));
cmin=min(SVP.ssp(:,2));
%
fprintf(fid1,'RAM PE RUN NAME = %s\n',ram_root);
fprintf(fid1,'%9.2f\t%9.2f\t%9.2f\t\t\t freq zs zr\n',[freq zs 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',0,hw);
fprintf(fid1,'-1  -1 \n');
% Print out profile:
ndepth=size(SVP.ssp,1);
PE_prof_write(fid1,SVP.ssp,ndepth,GAP,lambda,nlam_bot,nlam_false,freq);
fclose(fid1);
%
%%%
%
function PE_prof_write(fid1,svp_sec,ndepth,GAP,lambda,nlam_bot,nlam_false,freq)
%
% Writes out a geoacoustic profile for a particular range
%
% Print out depth, sound speed in water column:
fprintf(fid1,'%9.3f %9.3f \t\t\t\t\t z cw \n',svp_sec(1,1),svp_sec(1,2));
%	
for jz=2:ndepth,
	fprintf(fid1,'%9.3f %9.3f \n',svp_sec(jz,1),svp_sec(jz,2));
end
fprintf(fid1,'-1  -1 \n');
%
zbot=svp_sec(ndepth,1);
% Interpret negative sediment sound speed as ratio of csed/cw.  Bottom of first layer
% based on top value, so construct dc to add to cp2 if necessary (also add to subsequent
% layers since this might be a continuous layer, and dc probably small anyway):
% ([1500 m/s used in adi_svp_interpret as nominal value]);
csed=GAP.geo(1,1,1);
if csed<0, cs_cw_rat=-csed; csed=cs_cw_rat*svp_sec(ndepth,2); else, csed=1500; end
dc=csed-1500;
%
zdep=cumsum([0; GAP.h]);
% Print out depth, sound speed in bottom:
tabs='\t\t\t\t\t z cb ';
for jb=1:GAP.nlay,
	fprintf(fid1,['%9.3f %9.3f ',tabs,'\n'],zbot + zdep(jb),GAP.geo(1,1,jb) + dc); tabs=[];
	if jb<GAP.nlay, fprintf(fid1,'%9.3f %9.3f \n',zbot + zdep(jb+1),GAP.geo(2,1,jb) + dc); end
end
fprintf(fid1,'-1  -1 \n');
%
% Print out depth, density in bottom:
tabs='\t\t\t\t\t z rhob ';
for jb=1:GAP.nlay,
	fprintf(fid1,['%9.3f %9.3f ',tabs,'\n'],zbot + zdep(jb),GAP.geo(1,3,jb)); tabs=[];
	if jb<GAP.nlay, fprintf(fid1,'%9.3f %9.3f \n',zbot + zdep(jb+1),GAP.geo(2,3,jb)); end
end
fprintf(fid1,'-1  -1 \n');
%
% Print out depth, attenuation in bottom:
tabs='\t\t\t\t\t z attnb ';
% Compute frequency exponent factor so that fexp=1 results in ffac=1:
% Must also convert from dB/m/kHz to dB/lambda:
for jb=1:GAP.nlay,
	lamfac=.001*GAP.geo(1:2,1,jb);
	ffac=(freq/1000)^(GAP.fexp(1,jb)-1);	% disp(num2str([jb,GAP.fexp(1,jb),ffac])),
	attn_bot=lamfac(1)*ffac*GAP.geo(1,2,jb);
%     str=sprintf('jb = %i, orca attn = %3.3f, ram attn = %3.3f, c = %4.1f, ffac=%10.4g', ...
%         [jb,GAP.geo(1,2,jb),attn_bot,GAP.geo(1,1,jb),ffac]); disp(str);
	fprintf(fid1,['%9.3f %9.3f ',tabs,'\n'],zbot + zdep(jb),attn_bot);
	tabs=[];
	if jb<GAP.nlay, 
		fprintf(fid1,'%9.3f %9.3f \n',zbot + zdep(jb+1),lamfac(2)*ffac*GAP.geo(2,2,jb)); 
	end
end
% Print out depth, attenuation in absorbing bottom:
z_hsp=zbot + zdep(GAP.nlay+1) + nlam_bot*lambda;
z_false=z_hsp + nlam_false*lambda;
fprintf(fid1,'%9.3f %9.3f \t\t\t\t\t zb attnb false bottom \n',z_hsp,attn_bot);
fprintf(fid1,'%9.3f %9.3f \n',z_false,max(attn_bot,10));
fprintf(fid1,'-1  -1 \n');
%
%%%
%
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;

%
% 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.
