function[]=sphereEX2()
%Kevin Luna
%SpherEX1.m
%The computes the axis symetric r solution and derivative in spherical coordinates of the
%model problem (1/r^2)(\beta r^2 u_r)_r =f(r) + c\delta(x-\alpha) with 
%piecewise constant \beta using the IIM and performs a grid refinement analysis on them. 
%The structure of this code is the same as Ex1.m
% Change the functions in the 'Interface BVP functions' area to compute new
% solutions and derivatives. The ouputs are grid rifenment analysis plots
% of the solution and the derivative at the interface, and error plot for
% the last refinement, and approximate solution as a function of only r, and a
% 3D level surface plot at the interface in all three varibles

global alpha c betaminus betaplus W bf

n=20;
maxind=8;       % Number of refinements
hvec=zeros(maxind,1);       % Preallocation vector for h values
e1inf=zeros(maxind,1);      % Preallocation vector for error  values
T1=zeros(maxind,1);     % Preallocation vector for truncation error values
T2=zeros(maxind,1);     % Preallocation vector for truncation values






Ad=uDerexactminus(alpha);%exact left deriavative
B=uDerexactplus(alpha);%exact right  deriavative

%preallocation vectors for derivatives at alpha and errors there
UalphaprimeR=zeros(maxind,1);
UalphaprimeRE=zeros(maxind,1);
UalphaprimeL=zeros(maxind,1);
UalphaprimeLE=zeros(maxind,1);
Ualphaprime3=zeros(maxind,1);
Ualphaprime3E=zeros(maxind,1);
DUalphaprime3=zeros(maxind,1);
DUalphaprime3E=zeros(maxind,1);





for i=1:maxind
    [A,r,U,dr,F,gamma1,gamma2,Cj1,Cj2,j]=polarex(n); % calling IIM FD solver 
    u=zeros(n-1,1);% Preallocation vector for exact solution values
    for k=1:n-1
             u(k)=uexact(r(k));% filling u with exact solution values
    end
    
T1(i)=gamma1(1)*u(j-1)+gamma1(2)*u(j)+gamma1(3)*u(j+1)-F(j); % Truncation error at X(j)  for current number of grid divisions  
T2(i)=gamma2(1)*u(j)+gamma2(2)*u(j+1)+gamma2(3)*u(j+2)-F(j+1);% Truncation error at X(j+1) for current number of grid divisions
hvec(i)=dr ;% computing h for current number of grid points
e1inf(i)=norm(u-U,inf); %computing error for current number of grid points

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Computing correction terms for the derivatives at alpha
X=r;
B11d=[1 ,1, 1];
B21d=[X(j-1)-alpha , X(j)- alpha, (betaminus/betaplus)*(X(j+1)-alpha)];
B31d=[.5*(alpha-X(j-1))^2, .5*(alpha-X(j))^2 , .5*(betaminus/betaplus)*(X(j+1)-alpha)^2];
Bj1d=[B11d;B21d;B31d];
gammaRHS1d=[0;1;0];
gamma1d=Bj1d\gammaRHS1d; 

B12d=[1 ,1 ,1];
B22d=[(betaplus/betaminus)*(X(j)-alpha) , X(j+1)- alpha, (X(j+2)-alpha)];
B32d=[ .5*(betaplus/betaminus)*(X(j)-alpha)^2, .5*(X(j+1)-alpha)^2 , .5*(X(j+2)-alpha)^2 ];
Bj2d=[B12d;B22d;B32d];
gammaRHS2d=[0;1;0];
gamma2d=Bj2d\gammaRHS2d; 


Cj1d=gamma1d(3)*(W+(X(j+1)-alpha)*(c/betaplus)-.5*((X(j+1)-alpha)^2)*( -bf/betaplus + (2*c/(betaplus*alpha))    ));
Cj2d=gamma2d(1)*(-W-(X(j)-alpha)*(c/betaminus)+.5*((X(j)-alpha)^2)*( -(bf/betaminus) + (2*c/(betaminus*alpha))   ));






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Computing two-sided derivatives
%left
Ualphaprime3(i)=gamma1d(1)*U(j-1)+gamma1d(2)*U(j)+ gamma1d(3)*U(j+1)-Cj1d;
Ualphaprime3E(i)=abs(Ad-Ualphaprime3(i));

%right
DUalphaprime3(i)=gamma2d(1)*U(j)+gamma2d(2)*U(j+1)+ gamma2d(3)*U(j+2)-Cj2d;
DUalphaprime3E(i)=abs(B-DUalphaprime3(i));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Computing one-sided derivatives

UalphaprimeR(i)= -(3/(2*hvec(i)))*U(j+1) +(2/hvec(i))*U(j+2) -(1/(2*hvec(i)))*U(j+3);
UalphaprimeRE(i)=abs(B-UalphaprimeR(i));
UalphaprimeL(i)= (3/(2*hvec(i)))*U(j) -(2/hvec(i))*U(j-1) +(1/(2*hvec(i)))*U(j-2);
UalphaprimeLE(i)=abs(Ad-UalphaprimeL(i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%















n=2*n;
end


%%%OUTPUTS
% loglog plot of solution error

% figure(1)
% plot(log10(hvec),log10(e1inf),'or')
% legend('Immersed Interface Method','Location','northwest')
% xlabel('log_{10}(h)')
% ylabel('log_{10}( ||U-u||_\infty )')
% title('Solution Error Plot on a log-log Scale for the IIM')
% hold on
% p1=polyfit(log10(hvec),log10(e1inf),1); %fit line onto loglog data
% xh0=linspace(hvec(1),hvec(maxind),1000);
% x0=log10(xh0);
% y0=polyval(p1,x0);
% plot(x0,y0)
% axis([-6 -1 -6 -1])
% axis('square')
%             print('solnsph2','-depsc')
% hold off
% 
% % loglog plot of the left  derivative  error
figure(2)
plot(log10(hvec),log10(DUalphaprime3E),'or',log10(hvec),log10(UalphaprimeRE),'sq')
title('loglog Plot of U_x-(\alpha) Error')
xlabel('log_{10}(h)')
ylabel('log_{10}( ||U_x^-(\alpha)-u_x^-(\alpha)||_\infty )')
legend('Two-Sided Method(IIM)','One-sided method','Location','northwest')
hold on
p4=polyfit(log10(hvec),log10(DUalphaprime3E),1);
p6=polyfit(log10(hvec),log10(UalphaprimeRE),1);%fit line onto loglog data
xh1=linspace(hvec(1),hvec(maxind),1000);
x1=log10(xh1);
y1=polyval(p4,x1);
plot(x1,y1)
xh2=linspace(hvec(1),hvec(maxind),1000);
x2=log10(xh2);
y2=polyval(p6,x2);
plot(x2,y2)
M=log10([DUalphaprime3E(2:end);UalphaprimeRE(2:end)]);


axis([-6.5 -1 -6.5 -1])
axis('square')
        print('alphasph2','-depsc')
hold off
% 
% 
% %difference error plot
% figure(3)
% plot(X,U-u)
% title('Error (difference) Plot')
% xlabel('x' )
% ylabel('U-u(x)')
%         print('errorsph2','-depsc')



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p1=polyfit(log10(hvec),log10(e1inf),1); %fit line onto loglog data
fprintf('The slope of the  solution error plot on the  infinity norm plot is %.8f \n',p1(1))

p2=polyfit(log10(hvec),log10(T1),1); %fit line onto loglog data
fprintf('The slope of the truncation  error at j for the  loglog infinity norm data is %.8f \n',p2(1))


p3=polyfit(log10(hvec),log10(T2),1); %fit line onto loglog data
fprintf('The slope of the truncation  error at j+1 for the loglog infinity norm data is %.8f \n',p3(1))


p4=polyfit(log10(hvec),log10(Ualphaprime3E),1);
fprintf('The slope of the left two sided U''(alpha) error plot is %.8f \n',p4(1))


p5=polyfit(log10(hvec),log10(DUalphaprime3E),1);%fit line onto loglog data
fprintf('The slope of the right two sided U''(alpha) error plot is %.8f \n',p5(1))

p6=polyfit(log10(hvec),log10(UalphaprimeLE),1);%fit line onto loglog data
fprintf('The slope of the left  one sided U''(alpha) error plot is %.8f \n',p6(1))

p7=polyfit(log10(hvec),log10(UalphaprimeRE),1);%fit line onto loglog data
fprintf('The slope of the right  one sided U''(alpha) error plot is %.8f \n',p7(1))




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%local truncation error
n=n/2;
T=zeros(n-1,1);
T(1)= A(1,1)*u(1)+ A(1,2)*u(2) -F(1);
for i=2:n-2
    if i==j
     T(j)=gamma1(1)*u(j-1)+gamma1(2)*u(j)+gamma1(3)*u(j+1)-F(j); % Truncation error at X(j)  for current number of grid divisions  
    else if i==j+1
        T(j+1)=gamma2(1)*u(j)+gamma2(2)*u(j+1)+gamma2(3)*u(j+2)-F(j+1);
        end
    end
    
    
    
  T(i)=  A(i,i-1)*u(i-1) +A(i,i)*u(i) + A(i,i+1)*u(i+1)-F(i); 
end
T(n-1)=A(n-1,n-2)*u(n-2)+A(n-1,n-1)*u(n-1)-F(n-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solution plot
[M,i]=min(T)
% figure(4)
% plot(r,U,'o',r,u)
% title('Solution  and Approximation plot for \beta (r^2 u_r)_r = r^2 f(r) ')
% xlabel('r')
% ylabel('u(r)')
% legend('Approximation', 'Exact Solution')
%             print('plotsph2','-depsc')

%plotting 3-d contours
theta=linspace(0,pi,n-1) ;
phi=linspace(0, pi/2,n-1);
[theta,phi] = meshgrid(theta,phi);
 R1=U(j);
 R2=U(j+1);
 R3=U(j-n/4);
 R4=U(j+n/4);
 [x1,y1,z1]=sph2cart(theta,phi,R1);
 [x2,y2,z2]=sph2cart(theta,phi,R2);
 [x3,y3,z3]=sph2cart(theta,phi,R3);
 [x4,y4,z4]=sph2cart(theta,phi,R4);
 
%  figure(5)
%  k1=surf(x1,y1,z1);
%    axis([0 1 0 1 0 .75]) 
% hold on
% k2=surf(x2,y2,z2);
%    %axis([0 1  1 -1 1]) 
% set(k1,'edgecolor','none')
% set(k2,'edgecolor','none')
% title('Level Surface Before and After Interface')
% xlabel('x')
% ylabel('y')
% zlabel('z')


end






function[A,r,U,dr,F,gamma1,gamma2,Cj1,Cj2,j]=polarex(n)


global c alpha betaminus betaplus W bf
alpha=2/3;
betaminus=1;
betaplus=100;


c=ccalc(alpha); %[\beta u_x]_{x=\alpha} =c jump
bf=bfcalc(alpha); %[f] jump
W=bucalc(alpha);%  [u]=W jump

R2=1;
ua=uexact(0);
ub=uexact(R2);
[A,r,rhalf,dr]=spheregrid(n,R2);


X=r;
Y=zeros(n-1,1);
for i=1:n-1
Y(i)=irreg_j_index(X(i));
end
[M,j] = max(Y);


F=zeros(n-1,1); 
for i=1:n-1 
   
    F(i)=r(i)*r(i)*f(r(i)); 
    
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%IIM computation 


B11=[1 ,1, 1];
B21=[X(j-1)-alpha , X(j)- alpha, (betaminus/betaplus)*(X(j+1)-alpha)];
B31=[.5*(alpha-X(j-1))^2, .5*(alpha-X(j))^2 , .5*(betaminus/betaplus)*(X(j+1)-alpha)^2];
Bj1=[B11;B21;B31];
gammaRHS1=[0;2*alpha*betaminus;alpha*alpha*betaminus];
gamma1=Bj1\gammaRHS1; 

B12=[1 ,1 ,1];
B22=[(betaplus/betaminus)*(X(j)-alpha) , X(j+1)- alpha, (X(j+2)-alpha)];
B32=[ .5*(betaplus/betaminus)*(X(j)-alpha)^2, .5*(X(j+1)-alpha)^2 , .5*(X(j+2)-alpha)^2 ];
Bj2=[B12;B22;B32];
gammaRHS2=[0;2*betaplus*alpha;alpha*alpha*betaplus];
gamma2=Bj2\gammaRHS2; 


Cj1=gamma1(3)*(W+(X(j+1)-alpha)*(c/betaplus)-.5*((X(j+1)-alpha)^2)*( -bf/betaplus + (2*c/(betaplus*alpha))    ));
Cj2=gamma2(1)*(-W-(X(j)-alpha)*(c/betaminus)+.5*((X(j)-alpha)^2)*( -(bf/betaminus) + (2*c/(betaminus*alpha))   ));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Adding correction terms and changes to A
F(j)=F(j)+Cj1;
F(j+1)=F(j+1)+Cj2;

A(1:j-1,:)=betaminus.*A(1:j-1,:);
A(j+2:end,:)=betaplus.*A(j+2:end,:);


A(j,j-1)=gamma1(1);
A(j,j)=gamma1(2);
A(j,j+1)=gamma1(3);

A(j+1,j)=gamma2(1);
A(j+1,j+1)=gamma2(2);
A(j+1,j+2)=gamma2(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Adjusting boundary values based on discretization
F(1)=F(1)+ua/4 ;%%%%  ; 
F(n-1)=F(n-1) -(ub*betaplus*rhalf(end))/(dr*dr);

U=A\F ;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%INTERFACE BVP FUNCTIONS 
% Note that "-" means values BEFORE \alpha(interface point)
%AND that "+" means values AFTER \alpha(interface point)


function[y]=uexactminus(r)
%exact u- solution
global k1 k2
 y= cos(r^3);
end

function[y]=uexactplus(r)
% exact u+ solution
global k1 k2
 y = sin(r^2) ; %-.30;
end



function[y]=uDerexactminus(r)
%exact u_x-
y= -3*(r^2)*sin(r^3);
 
end


function[y]=uDerexactplus(r)
%exact u_x+

 y= 2*r*cos(r^2);
end


function[y]=fminus(r)
%f- RHS of ODE
global betaminus
y= -3*r*(4*sin(r^3) + 3*(r^3)*cos(r^3));
y=betaminus*y;
end

function[y]=fplus(r)
%f+ RHS of ODE
global betaplus
 y= 6*cos(r^2)- 4*(r^2)*sin(r^2);
y=betaplus*y;
end

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

function[RHS]=f(x)
% f(x)function
global alpha 
if x <= alpha
    RHS=fminus(x); 
    
else 
        
   RHS=fplus(x);
   
end
end

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

function [y] = irreg_j_index( x )
%outputs values of x <= \alpha and sets values of x > \alpha to 0
global alpha
if x <= alpha
    y=x;  
else   
  y=0 ;  
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[v]=ccalc(x) 
global betaplus betaminus 
%[\beta u_x] computation
v=betaplus*uDerexactplus(x)-betaminus*uDerexactminus(x);  
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=bfcalc(x)
%[f] computation
y= fplus(x) - fminus(x) ; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=uexact(x)
%Exact solution function
global alpha
if x <= alpha
    y=uexactminus(x); 
else 
    y=uexactplus(x); 
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[y]=bucalc(x)
%[u] computation
y=uexactplus(x)-uexactminus(x);
end