function ls52(action)
%written by John Sylvester, Mathematics Department,U. of Washington
%sylvest@math.washington.edu
%last edited 7/10/98
% 9/28/98 fixed problem in a2n -- generated nan's when c = constant

 global hr ylims low high
 global r rr   w weight absweight stopp ls51fig

 fname = mfilename;

if nargin<1
	action = 'initialize';
end
%====================================
%====================================
if strcmp(action,'initialize')	
%====================================
%====================================
	low = 0;
	high = inf;
    ls51fig =figure( ...
        'Name','LayerStrip', ...
		'backingstore','off',...
        'NumberTitle','off',...
		'Visible','off',...
		'interruptible','on');
    axrndl=axes( ...
        'Units','normalized', ...
        'tag','raxis',...
        'Position',[0.05 0.05 0.60 0.40]);
	title('Modulus of Reflection Coeficient');
	zoom on
	axalphahndl=axes( ...
        'Units','normalized', ...
        'tag','alphaaxis',...
        'Position',[0.05 0.55 0.60 0.40]);
	title('Alpha')
	zoom on
%===================================
    % Information for all buttons
%====================================

    top=0.95;
    bottom=0.05;
    labelColor=[0.8 0.8 0.8];
    yInitPos=0.90;
    left=0.70;
    btnWid=0.25;
    btnHt=0.05;
% Spacing between the button and the next command's label
    spacing=0.02;
    
%====================================
    % The CONSOLE frame
%====================================

    frmBorder=0.02;
    frmPos=[left-frmBorder bottom-frmBorder btnWid+3*frmBorder 0.9+2*frmBorder];
    h=uicontrol( ...
        'Style','frame', ...
        'Units','normalized', ...
        'Position',frmPos, ...
	'BackgroundColor',[0.5 0.5 0.5]);
	
%====================================
    % The Input R button
%====================================
    btnNumber=1;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Load R';
    callbackStr=[fname '(''inputr'')'];
    
    % Generic button information
    btnPos=[left-.01 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
%====================================
    % The Output R button
%====================================
    btnNumber=1;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Save R';
    callbackStr=[fname '(''outputr'')'];
    
    % Generic button information
    btnPos=[left+btnWid/4 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);	
%====================================
    % The Input Alpha button
%====================================
    btnNumber=1;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Load A';
    callbackStr=[fname '(''inputalpha'')'];
    
    % Generic button information
    btnPos=[left+btnWid/2+.01 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
%====================================
    % The Output  Alpha button
%====================================
    btnNumber=1;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Save A';
    callbackStr=[fname '(''outputalpha'')'];
    
    % Generic button information
    btnPos=[left+3*btnWid/4+.02 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);	

%====================================
    % The Stop button
%====================================
    btnNumber=1.9;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Stop';
    callbackStr= [fname '(''stop'')'];
    
    % Generic button
    btnPos=[left-.01 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
        
%====================================
    % The Grid button
%====================================
    btnNumber=1.9;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Grid';
    callbackStr=[fname '(''grid'')'];

  % Generic button
    btnPos=[left+ btnWid/4 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
%====================================
    % The Hold Alpha button
%====================================
    btnNumber=1.9;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Hold A';
    callbackStr='';

  % Generic button
    btnPos=[left+ btnWid/2+.01 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','check', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'tag','holdA',...
        'Callback',callbackStr);

%====================================
    % The Hold R button
%====================================
    btnNumber=1.9;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Hold R';
    callbackStr='';

  % Generic button
    btnPos=[left+ 3*btnWid/4+.02 yPos-spacing btnWid/4 btnHt];
    uicontrol( ...
        'Style','check', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'tag','holdr',...
        'Callback',callbackStr);

 
 %====================================  
 %%% The Genarate Reflection Coefficient Button
%====================================
    btnNumber=2.8;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['Alpha2R'];
    callbackStr = [fname '(''forward'')'];

    % Generic button information
    btnPos=[left-.01 yPos-spacing btnWid/3 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
%====================================
    % The Layerstrip button
%====================================
    btnNumber=2.8;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='R2Alpha';
    callbackStr=[fname '(''layerstrip'')'];
    
    % Generic button information
    btnPos=[left+btnWid/3+.01 yPos-spacing btnWid/3 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
 %====================================  
 %%% The alpha2n button
%====================================
    btnNumber=2.8;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['a2n'];
    callbackStr = [fname '(''a2n'')'];

    % Generic button information
    btnPos=[left+2*btnWid/3+.02 yPos-spacing btnWid/3 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Callback',callbackStr);
%====================================
    % The Energy Label
%====================================

    btnNumber=3.7;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['Energy0' blanks(10) 'Energy'] ;
    callbackStr = '';

    % Generic button
    btnPos=[left yPos btnWid btnHt/2];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor',[0.5 0.5 0.5], ...
        'ForegroundColor','white', ...
        'String',labelStr, ...
        'Callback',callbackStr);
		
%====================================
 % The Energy0 Text
    btnNumber=4.7;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=num2str(0);
    callbackStr = '';

    % Generic button 
    btnPos=[left yPos+spacing btnWid/2 btnHt];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','energy0',...
        'Callback',callbackStr);
%====================================
%====================================
 % The Energy Text
    btnNumber=4.7;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=num2str(0);
    
    % Generic button 
    btnPos=[left+btnWid/2+.01 yPos+spacing btnWid/2 btnHt];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','energy',...
        'Callback',callbackStr);
%====================================
%====================================
    % The Stepsize and Depth Label
    btnNumber=5.6;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['Stepsize' blanks(20) 'Depth'];
    callbackStr = '';

    % Generic button 
    btnPos=[left yPos+spacing btnWid btnHt];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor',[0.5 0.5 0.5], ...
        'ForegroundColor','white', ...
        'String',labelStr, ...
        'Callback',callbackStr);
		
%====================================
 % The stepsize edit
    btnNumber=5.5;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
	
    labelStr= '0.5';
    callbackStr = '';
    
    % Generic button information
    btnPos=[left yPos-spacing btnWid/2 btnHt];
    uicontrol( ...
        'Style','edit', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','stepsize',...
        'Callback',callbackStr);
%====================================
 % The depth edit
    btnNumber=5.5;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
	
    labelStr= '40.0';
    callbackStr = '';
    
    % Generic button
    btnPos=[left+ btnWid/2+.01  yPos-spacing btnWid/2 btnHt];
    depthHndl=uicontrol( ...
        'Style','edit', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','depth',...
        'Callback',callbackStr);
%====================================
% The number of DFT pts and Scale title
%====================================
    btnNumber=6.5;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['DFTpoints', blanks(15), 'w Scale Factor'];
    callbackStr = '';

    % Generic button 
    btnPos=[left yPos-spacing btnWid btnHt];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor',[0.5 0.5 0.5], ...
        'ForegroundColor','white', ...
        'String',labelStr, ...
        'Callback',callbackStr);
		
%====================================
 % The DFT  edit
%====================================
    btnNumber=7.3;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
	
    labelStr= '2^8';
    callbackStr = '';
    
    % Generic button
    btnPos=[left yPos btnWid/2 btnHt];
    uicontrol( ...
        'Style','edit', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','dft',...
        'Callback',callbackStr);
%====================================
 % The Scale Edit
%====================================
    btnNumber=7.3;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
	
    labelStr= '1';
    callbackStr = '';
    
    % Generic button
    btnPos=[left+ btnWid/2+.01  yPos btnWid/2 btnHt];
   uicontrol( ...
        'Style','edit', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','scale',...
        'Callback',callbackStr);
%====================================
    % The reflection coefficient Formula Label
%====================================
    btnNumber=8.3;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['Input Formula for R'] ;
    callbackStr = '';

    % Generic button
    btnPos=[left-.01 yPos btnWid+.01 btnHt];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor',[0.5 0.5 0.5], ...
        'ForegroundColor','white', ...
        'String',labelStr, ...
        'Callback',callbackStr);
		
%====================================
 % % The Reflection Coefficient Formula
%====================================
    btnNumber=8.5;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
	
    labelStr= '.7*exp(2*i*w*5)./(1-i*w).^2';
    callbackStr = [fname '(''rformula'')'];
    
    % Generic button information
    btnPos=[left-.01 yPos-spacing btnWid+.04 btnHt];
	uicontrol( ...
        'Style','edit', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','rformula',...
        'Callback',callbackStr);
        
%====================================
    % The alpha Formula Label
%====================================
    btnNumber=9.4;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr=['Input Formula for Alpha'] ;
    callbackStr = '';

    % Generic button
    btnPos=[left-.01 yPos btnWid+.01 btnHt/2];
    uicontrol( ...
        'Style','text', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor',[0.5 0.5 0.5], ...
        'ForegroundColor','white', ...
        'String',labelStr, ...
        'Callback',callbackStr);
		
%====================================
 % % The Alpha Formula
%====================================
    btnNumber=9.9;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
	
    labelStr= '.5*((x>6).*(x<8)-(x>10).*(x<12))';
    callbackStr = [fname '(''alphaformula'')'];
    
    % Generic button information
    btnPos=[left-.01 yPos-spacing btnWid+.04 btnHt];
    uicontrol( ...
        'Style','edit', ...
        'Units','normalized', ...
        'Position',btnPos, ...
		'BackgroundColor','white', ...
        'ForegroundColor','black', ...
        'String',labelStr, ...
        'tag','Aformula',...
        'Callback',callbackStr);


%====================================
    % ThePayley-Wiener Popoup 
%====================================
    btnNumber= 11;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='R | R/(1-abs(R)^2)| 1 /(1-abs(R)^2)-1';
    callbackStr=[fname '(''paleywiener'')'];
    
    % Generic button information
    btnPos=[left yPos-spacing btnWid/2 btnHt];
    uicontrol( ...
        'Style','popupmenu', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'tag','PW',...
        'Callback',callbackStr);
%====================================
    % The Linear/Nonlinear  Popup
%====================================
    btnNumber=11;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='Nonlinear | Linear';
    callbackStr='';
    
    % Generic button information
    btnPos=[left+ btnWid/2+.01 yPos-spacing btnWid/2 btnHt];
    uicontrol( ...
        'Style','popupmenu', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'tag','linear',...
        'Callback',callbackStr);
%====================================
    % The R- plot Popup
%====================================
    btnNumber=12;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr='plot abs(r) | real and imag R ';
    callbackStr=[fname '(''plotr'')'];
    
    % Generic button information
    btnPos=[left yPos-spacing btnWid/2 btnHt];

    uicontrol( ...
        'Style','popupmenu', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'tag','plotr',...
        'Callback',callbackStr);
		
%====================================
    % The Bandwidth limiting button
%====================================
    btnNumber=12;
    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
    labelStr= 'Bandwidth';
    callbackStr=[fname '(''bandwidth'')'];
    
%    % Generic button information
    btnPos=[left+btnWid/2+.01 yPos-spacing btnWid/2 btnHt];
    uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'tag','smoothing', ...
        'Callback',callbackStr);

%====================================
    % The Smoothing  Popup
%====================================
%    btnNumber=12;
%    yPos=top-btnHt-(btnNumber-1)*(btnHt+spacing);
%    labelStr='Smooth off| Smooth on';
%    callbackStr='';
%    
%    % Generic button information
%    btnPos=[left+btnWid/2+.01 yPos-spacing btnWid/2 btnHt];
%    uicontrol( ...
%        'Style','popupmenu', ...
%        'Units','normalized', ...
%        'Position',btnPos, ...
%       'String',labelStr, ...
%        'tag','smoothing', ...
%       'Callback',callbackStr);
%		
%===========================================
    % The CLOSE button
%====================================
    labelStr='Close';
    callbackStr=[fname '(''close'')'];
    uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'Position',[left bottom btnWid btnHt], ...
        'String',labelStr, ...
        'Callback',callbackStr);


 
 %====================================   
    figure(ls51fig);

   


%====================================
%====================================
elseif strcmp(action,'layerstrip')
%====================================
%====================================
	axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	axrhndl = findobj(ls51fig,'tag','raxis');
	energyHndl= findobj(ls51fig,'tag','energy');
	energy0Hndl = findobj(ls51fig,'tag','energy0'); 
	deltaHndl=findobj(ls51fig,'tag','stepsize'); 
	depthHndl= findobj(ls51fig,'tag','depth');
	smoothHndl = findobj(ls51fig,'tag','smoothing'); 
	linearHndl = findobj(ls51fig,'tag','linear'); 
	plotrHndl =findobj(ls51fig,'tag','plotr'); 



depth = str2num(get(depthHndl,'String'));
delta = str2num(get(deltaHndl,'String'));
smoothing = get(smoothHndl,'value');
linear = get(linearHndl,'Value')-1;

if linear == 1
	nitermax = 1;
else
	nitermax = 25;
end

%initialize
Nalpha = floor(depth/delta) +1;
alpha = zeros(1,Nalpha);

%compute shifts and Fourier transform of the Heavyside
twoiw = 2*i*w;
s = exp(twoiw*delta);
sminus = 1./s;
Heavyhat = (1-sminus)./twoiw/2;


eval([fname '(''prepalphaplot'')'])
ylims = get(axalphahndl,'Ylim');


		x = (0:Nalpha)*delta;
		halpha = line(x(1),0,'erasemode','background','tag','alphaplot');
		set(axalphahndl,'Xlim',[0,depth],'drawmode','fast');
		title('Alpha')
		  


factortol = 1e-5;
[R,A,niter] = factor(r,weight,absweight,factortol,nitermax,smoothing);

eR0 = energy(R,absweight);
set(   energy0Hndl,'String',num2str(  eR0  )   )
drawnow

n = 0;
eR= eR0;

%etol = 1e-3;
stopp = 0;
   while n<Nalpha  & stopp == 0  %%& eR/eR0>etol
	n = n+1;
	smoothing = get(smoothHndl,'value');
	[R,A,niter] = ...
		factor(R.*sminus,weight,absweight,factortol,nitermax,smoothing);
	eR = energy(R,absweight);
	alpha(n) = A/Heavyhat;
	%%plot stuff
		set(halpha,'Xdata',x(1:n),'Ydata',alpha(1:n))
		if alpha(n) < ylims(1) 
				ylims(1) = 2*alpha(n);
				set(axalphahndl,'Ylim',ylims)
		elseif alpha(n) > ylims(2)
				 ylims(2) = 2*alpha(n);
				 set(axalphahndl,'Ylim',ylims)
		end
		if get(plotrHndl,'Value') == 1
			set(hr,'Ydata',abs(R))
		elseif get(plotrHndl,'Value') == 2
			set(hr(1),'Ydata',real(R))
			set(hr(2),'Ydata',imag(R))
		end
		set(  energyHndl,'String',num2str(  eR  )   )
		drawnow

	%%end plot stuff
   end
   set(halpha,'Xdata',x,'Ydata',[alpha,zeros(1,length(x)-length(alpha))])
   
%====================================
%====================================
elseif strcmp(action,'rformula')
%====================================
%====================================
     rformHndl= findobj(ls51fig,'tag','rformula');
     scaleHndl  = findobj(ls51fig,'tag','scale');
     DFTHndl = findobj(ls51fig,'tag','dft');

    N = str2num(get(DFTHndl,'String'));
    scalestr = get(scaleHndl,'string');
    [w,weight,absweight] = makefrequencies(N);
    rformula = strrep(get(rformHndl,'String'),'w',['(w/' scalestr ')']);
    eval(['rr =' rformula ';'])
    eval([fname '(''paleywiener'')'])
    
%====================================
elseif strcmp(action,'alphaformula')
%====================================
%====================================
	 axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	 deltaHndl=findobj(ls51fig,'tag','stepsize');
	 depthHndl= findobj(ls51fig,'tag','depth');
	 alphaHndl = findobj(ls51fig,'tag','Aformula');

 	delta = str2num(get(deltaHndl,'String'));
    depth = str2num(get(depthHndl,'String'));
    x = 0:delta:depth;
    eval(['alpha =' get(alphaHndl,'String') ';'])
    
	eval([fname '(''prepalphaplot'')'])
	
 	halpha = line(x,alpha,'color','b','tag','alphaplot');

%====================================
%====================================
elseif strcmp(action,'inputr')
%====================================
%====================================
     DFTHndl = findobj(ls51fig,'tag','dft');
     scaleHndl  = findobj(ls51fig,'tag','scale');

     [name, path] = uigetfile('*.*at');
     if name ~= 0
       currentdir = pwd;
       cd(path);
       load(name);
       num = findstr(name,'.')-1;
       basename = name(1:num);
	   eval(['NN = size(', basename,');'])
	   if min(NN)>1
	   		%%need to interpolate to w grid
			if min(NN) == 3
				%% if real and imaginary parts are in 
				%% columns two and three, put the sum in the 
				%% second column and get rid of the third
				eval(['wdata =',basename,'(:,1);'])
				eval(['rdata =',basename,'(:,2)+i*',basename,'(:,3);'])
			else
				eval(['wdata =',basename,'(:,1);'])
				eval(['rdata =',basename,'(:,2);'])
			end
			
			if min(wdata)>0
				rdata = [ flipud(conj(rdata));rdata ];
				wdata = [ - flipud(wdata); wdata ];
			elseif min(wdata)==0
				l=length(rdata);
				rdata = [ flipud(conj(rdata(2:l)));rdata ];
				wdata = [ - flipud(wdata(2:l)); wdata ];
			end
			
			%%Now scale w's for better interpolation
			
			wdata = str2num(get(scaleHndl,'string'))*wdata;
			
	   		N = str2num(get(DFTHndl,'String'));
            		[w,weight,absweight] = makefrequencies(N);
			ind = find(w<max(wdata) & w>min(wdata));
			rr = zeros(size(w));
			rr(ind) = interp1(wdata,rdata,w(ind));		
	   else	
          	eval(['rr =' basename ';']);
		rr = rr(:).';
        	N = length(rr)/2;
		set(DFTHndl,'String',num2str(N))
        	[w,weight,absweight] = makefrequencies(N);
	   end
        cd(currentdir)
        eval([fname,'(''paleywiener'')'])
    end
    
%====================================
%====================================
elseif strcmp(action,'inputalpha')
%====================================
%====================================
	 axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	 deltaHndl=findobj(ls51fig,'tag','stepsize');

     [name, path] = uigetfile('*.mat');
     if name ~= 0
        currentdir = pwd;
        cd(path);
        load(name);
        cd(currentdir)
        num = findstr(name,'.mat')-1;
        basename = name(1:num); 
        eval(['alpha =' basename ';']);
        Nalpha =  length(alpha);
        delta = str2num(get(deltaHndl,'string'));
        x = delta*(0:Nalpha-1);
	
        eval([fname '(''prepalphaplot'')'])
 
	halpha = line(x,alpha,'tag','alphaplot');
	 end
    
%====================================
%====================================
elseif strcmp(action,'outputr')
%====================================
%====================================
	 %axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	 %halpha = get(axalphahndl,'Userdata');  

     [name, path] = uiputfile('*.mat','Save As');
    if name ~= 0
       currentdir = pwd;
       cd(path);
       num = findstr(name,'.mat')-1;
       basename = name(1:num);
       %alphaname = ['alpha', basename];
       %eval([alphaname, '= get(halpha,''Ydata'');']);
       eval([basename,'= rr;']);
       %save(name,basename,alphaname);
	   save(name,basename);
       cd(currentdir)
    end
    
%====================================
%====================================
elseif strcmp(action,'outputalpha')
%====================================
%====================================
	 axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	 halpha =  findobj(ls51fig,'tag','alphaplot');

     [name, path] = uiputfile('*.mat','Save As');
   if name ~= 0     
     currentdir = pwd;
     cd(path);
     num = findstr(name,'.mat')-1;
     basename = name(1:num);
     eval([basename, '= get(halpha,''Ydata'');']);
     save(name,basename);
    
     cd(currentdir)
   end
%====================================
%====================================
elseif strcmp(action,'forward')
%====================================
%====================================
     DFTHndl = findobj(ls51fig,'tag','dft');
	 alphaHndl = findobj(ls51fig,'tag','Aformula');
	 axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	 halpha =   findobj(ls51fig,'tag','alphaplot');
	 deltaHndl=findobj(ls51fig,'tag','stepsize');
	 depthHndl= findobj(ls51fig,'tag','depth');
     linearHndl = findobj(ls51fig,'tag','linear'); 

    N = str2num(get(DFTHndl,'String'));
    [w,weight,absweight] = makefrequencies(N);


    % compute forward map


 % if strcmp(get(alphaHndl,'String'),'screen')
    x = get(halpha, 'Xdata');
    delta = x(2) -x(1);
    alpha = get(halpha, 'Ydata');

 % else
 %   delta = str2num(get(deltaHndl,'String'));
 %   depth = str2num(get(depthHndl,'String'));
 %    x = 0:delta:depth;
 %   eval(['alpha =' get(alphaHndl,'String') ';'])
 % end
     linear = get(linearHndl,'Value')-1;
     rr = rf(alpha,w,delta,linear);
     eval([fname,'(''paleywiener'')'])

%====================================
%====================================
elseif strcmp(action,'plotr')
%====================================
%====================================
	 axrhndl = findobj(ls51fig,'tag','raxis');
     plotrHndl =findobj(ls51fig,'tag','plotr'); 
     energy0Hndl = findobj(ls51fig,'tag','energy0'); 

        axes(axrhndl)
        if get(findobj(ls51fig,'tag','holdr'),'value')==1
			if exist('hr')==1,set(hr,'linestyle',':'),end
		else
			cla
		end
		if get(plotrHndl,'Value') == 1
			hr =line(w,abs(r));
			title('Abs(R)')
		elseif get(plotrHndl,'Value') == 2
			hr =[line(w,real(r),'color','k'),line(w,imag(r),'color','r')];
			title('Real and Imag R')
		end
		set(axrhndl,'drawmode','fast','Xlim',[-5,5],'YlimMode','auto')
		set(hr,'erasemode','background');
		set( energy0Hndl,'String',num2str(  energy(r,absweight) ) )
		
		

%====================================
%====================================
elseif strcmp(action,'stop')
%====================================
%====================================

        stopp = 1;

%====================================
%====================================
elseif strcmp(action,'paleywiener')
%====================================
%====================================
    pwHndl = findobj(ls51fig,'tag','PW'); 

	r = rr.*(abs(w)>=low&abs(w)<=high);
	if get(pwHndl,'Value') == 2
		r = r./(1-abs(r).^2);
	elseif get(pwHndl,'Value') == 3
		r = (abs(r).^2) ./(1-abs(r).^2);
	end
	eval([fname,'(''plotr'')'])

%====================================
%====================================
elseif strcmp(action,'grid')
%====================================
%====================================
axalphahndl = findobj(ls51fig,'tag','alphaaxis');
axrhndl = findobj(ls51fig,'tag','raxis');

	if strcmp(get(axalphahndl,'Xgrid'),'on')
		set(axalphahndl,'Xgrid','off','Ygrid','off')
		set(axrhndl,'Xgrid','off','Ygrid','off')
	else
		set(axalphahndl,'Xgrid','on','Ygrid','on')
		set(axrhndl,'Xgrid','on','Ygrid','on')
	end

%====================================
%====================================
elseif strcmp(action,'close')
%====================================
%====================================
	
clear global r rr hr w weight absweight stopp
close(gcf)
%====================================
%====================================
elseif strcmp(action,'a2n')
%====================================
%====================================
	 axalphahndl = findobj(ls51fig,'tag','alphaaxis');
	 halpha = findobj(ls51fig,'tag','alphaplot');
	 deltaHndl = findobj(ls51fig,'tag','stepsize'); 
	 delta = str2num(get(deltaHndl,'String'));
	 alpha = get(halpha,'Ydata');
	
c = [1, cumprod(exp(-alpha .* delta))];
deltac = diff(c);

%if alpha = 0, c is constant so
deltaz = delta.* c(2:length(c));

%if alpha ~= 0, 
index = find(alpha~=0);
%c is an exponential with exponent alpha
deltaz(index) = -deltac(index)./alpha(index);

z =  [0,cumsum(deltaz)];

n = 1./c;
 %    c = [1, cumprod(exp(-alpha * delta))];
 %    deltaz = -diff(c)./alpha;
 %    z =  [0,cumsum(deltaz)];
 %    n = 1./c;

	figure, plot(z,n)
	title('Index of Refraction vs Physical Depth')
	xlabel('z'),ylabel('n')
	grid on, zoom on
	
%====================================
%====================================
elseif strcmp(action,'prepalphaplot')
%====================================
%====================================
	 	axalphahndl = findobj(ls51fig,'tag','alphaaxis');
		axes(axalphahndl)
		if get(findobj(ls51fig,'tag','holdA'),'value')==1
			flag = findobj(ls51fig,'tag','alphaplot');
			if ~isempty(flag)
				set(flag,'color','r','tag','')
			end
		else
			cla
			set(axalphahndl,'YlimMode','auto')
		end
%====================================
%====================================
elseif strcmp(action,'bandwidth')
%====================================
%====================================
		zoom off
		[vx, vy ,butt] = ginput(2);
		ss = length(vx);
		if ss == 0
			low = 0;
			high = inf;
		elseif ss == 1
			if butt == 1
				low = 0 ;
				high = vx(1);
			else

				low = vx(1);
				high = inf;
			end
		else
			low = vx(1,1);
			high = vx(2,1);
		end
		zoom on
		eval([fname,'(''paleywiener'')'])
end

%====================================
%=======makefrequencies.m
%====================================

function [w,weight,absweight] = makefrequencies(N);
thetaover2= ((-N:N-1) + 0.5)*pi/N/2;
w = tan(thetaover2);
weight =  i*w-1;
absweight = 1 + w.*conj(w);

%====================================
%=======rf.m
%====================================
function r = rf(alpha,w,delta,linear)
%compute  forward  scattering transform
%for piecewise constant alpha
%constant on intervals of length delta 
%compute at frequencies w

if nargin < 4 , linear = 0; end
Nalpha = length(alpha);
Nw = length(w);
r = zeros(1,Nw);
twoiw = 2*i*w;
s = exp(twoiw*delta);
sminus = 1./s;
Heavyhat = (1-sminus)./twoiw/2;
Heavyhat(w==0)=delta*ones(1,sum(w==0));

if linear == 0
	for n = Nalpha:-1:1
		r= s .* star(r,alpha(n)*Heavyhat);
	end
else
	for n = Nalpha:-1:1
		r= s .* (r + alpha(n)*Heavyhat);
	end
end
%====================================
%=======factor.m
%====================================
function [r,A,niter,energyerror] = factor(f,weight,absweight,tol,nitermax,smoothing)
%factor F = star(r,A)

if nargin < 6, smoothing=1; end
if nargin < 5, nitermax = 20; end
if nargin < 4, tol = 1e-3; end
N = length(f)/2;
r = zeros(size(f));
A = zeros(size(f));
nonlin = nitermax ~= 1;


niter = 0;
residual = 1;
while abs(residual)>=tol & niter<nitermax
	niter = niter+1;
	F = f.*(1+ conj(A).*r);
	c = fftshift(fft(fftshift(F.*weight)));
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%C = [zeros(1,N) , c(N+1:2*N)];
	%r = fftshift(ifft(fftshift(C)))./weight;
	%A = F-r;
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	C = [c(1:N),zeros(1,N)];
	A = fftshift(ifft(fftshift(C)))./weight;
	r = F-A;
	residual = norm(f- star(r,A),2)/norm(f,2);
end 
if smoothing == 2
r(2:2*N-1) = (2*r(2:2*N-1) +  r(1:2*N-2)+  r(3:2*N))/4;
end
if niter == nitermax & nonlin
	 disp(['Maximum number of iterations reached ',num2str(niter)])
end
%e = energy(f,absweight);
%energyerror = (e-energy(r,absweight)-energy(A,absweight))/e;
%====================================
%=======energy.m
%====================================
function e= energy(f,absweight)


	e =  -log(1-abs(f).^2)*(absweight.')/length(absweight);

%====================================
%=======star.m
%====================================
function F = star(r,b)
%conformal map of disk to disk mapping 0 to b
%F = (r+b)./(1+ conj(b).*r);

F = (r+b)./(1+ conj(b).*r);
%====================================
%=======alpha2n.m
%====================================
function [n,z] = alpha2n(alpha,x)

%output  is a list of values of n at poitns z; (1/n is piecewise linear)
%input is a list of values of piecewise constant alpha for x_j<x<x_j+1
%

deltax = diff(x);


c = [1, cumprod(exp(-alpha .* deltax))];
deltac = diff(c);

%if alpha = 0, c is constant so
deltaz = deltax.* c(2:length(c));

%if alpha ~= 0, 
index = find(alpha~=0);
%c is an exponential with exponent alpha
deltaz(index) = -deltac(index)./alpha(index);

z =  [x(1),cumsum(deltaz) + x(1)];

n = 1./c;

%====================================
%=======n2alpha.m
%====================================

function [Alpha,X] = n2alpha(n,z,dx)

%input  is a list of values of n at poitns z; (1/n is piecewise linear)
%output is a list of values of piecewise constant alpha for x_j<x<x_j+1
%


c = n(1) ./ n;   %speed of light is 1 on top---this is choice of units

%must remove deltac ==0's
%index = [1,find(diff(c)~=0)+1];
%c = c(index);
%z = z(index);

deltaz = diff(z);
deltac = diff(c);

N = length(c);

alpha = - deltac./deltaz;

%if alpha = 0, c is constant so
%deltax = deltaz./ c(2:length(deltax));

%if alpha ~= 0, 
index = find(alpha~=0);
Index = find(alpha ==0);

%c is an exponential with exponent alpha
deltax(Index) = deltaz(Index)./c(Index+1);
deltax(index) = log(  c(index)./c(index+1)  )./alpha(index);
x = [z(1),cumsum(deltax) + z(1)];

%Now interpolate to equispaced travel time (dx)
X = 0:dx:x(length(x));
Alpha = zeros(size(X));
for i = 1:length(x)-1
	I = find(X>= x(i) & X<x(i+1));
	Alpha(I) = alpha(i)*ones(size(I));
end