[OP,rhs]=bc_1d(M,bc_type,bc_vals,X)

function [OP,rhs]=bc_1d(M,bc_type,bc_vals,X)

% Usage: [OP,rhs]=bc_1d(integer,
% {bc_type1,...,bc_typeM},[bc_val1,...,bc_valM,[Left endpoint, Right
% endpoint])

% Returns 
%  OP: Matrix with Tau lines
%  rhs: vector with boundary conditions

% Inputs:
% M - discretization size of the Matrix
% bc_type -> See below 
% bc_vals -> bcs to be enforced
% X -> left and right interval

% Generate the tau line boundary conditions and the appropriate rhs vector

% Plus/minus conditions are no longer supported

% Inputs
% bc_type : vector of strings containing the bc type
%   Dirichlet
%   D_R = Diriclet boundary condition on the right
%   D_L = Diriclet boundary condition on the left
%   D_P = Summation of dirichlet conditions
%   D_M = Difference of dirichlet conditions
%   Neumann
%   N_R = Neumann on the right
%   N_L = Neumann on the left
%   N_P = Summation of neumann conditions
%   N_M = Difference of neumann conditions
%   Stress Free
%   SF1_R = Stress free for stream functions on the right
%   SF1_L = Stress free for stream functions on the left
%   SF2_R = Stress free for azimuthal velocity on the right
%   SF2_L = Stress free for azimuthal velocity on the left
%   SF3_R = Stress free for potential functions on the right
%   SF3_L = Stress free for potential functions on the left
%   D3_R = Third Derivative on the right
%   D3_L = Third Derivative on the left
%   D3_P = Summation of Third Derivative conditions
%   D3_M = Difference of Third Derivative conditions
%   I_RRR = Integral constraint on Third Derivative
%   I_DZ = Integral constraint on First Derivative in Z
%   I_SF2_tilde = Integral Constraint for SF2 condition


% Check that the appropriate conditions were passed
if (length(bc_type) ~= length(bc_vals))
   disp(sprintf('   Error in bc_1d.m : # of bc_types and bc_vals do not match'))
   return
end

len=length(bc_type);

RO=X(2);
RI=X(1);
F1=.5*(RO-RI);
G1=.5*(RO+RI);
    
OP=zeros(M);
rhs=zeros(M,1);


for l=1:len
   typ=bc_type(l);
   bc_pos=M-(len-l);
   if strcmp(typ,'D_R')
       for j=1:M
           OP(bc_pos,j)=1;
       end
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'D_L')
       for j=1:M
           OP(bc_pos,j)=1*(-1)^(j+1);
       end
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'N_R')
       for j=1:M
           OP(bc_pos,j)=(j-1)^2;
       end
       rhs(bc_pos)=bc_vals(l)*F1;
   elseif strcmp(typ,'N_L')
       for j=1:M
           OP(bc_pos,j)=(j-1)^2*(-1)^(j);
       end
       rhs(bc_pos)=bc_vals(l)*F1;
   elseif strcmp(typ,'D2_R')  % second derivative
       for j=1:M
           OP(bc_pos,j)=(1/3)*(j-1)^2*((j-1)^2-1);
       end
       rhs(bc_pos)=bc_vals(l)*F1^2;
   elseif strcmp(typ,'D2_L')
       for j=1:M
           OP(bc_pos,j)=(1/3)*((j-1)^2*((j-1)^2-1))*(-1)^(j+1);
       end
       rhs(bc_pos)=bc_vals(l)*F1^2;  
   elseif strcmp(typ,'SF1_R')  %stress free for stream function at right psi_rr - (psi_r)/r
       for j=1:M
           OP(bc_pos,j)=(1/3)*(j-1)^2*((j-1)^2-1)-(F1/RO)*((j-1)^2);  %RO
       end
       rhs(bc_pos)=bc_vals(l)*F1^2;
   elseif strcmp(typ,'SF1_L') %stress free for stream function at left
       for j=1:M
           OP(bc_pos,j)=((1/3)*(j-1)^2*((j-1)^2-1))*(-1)^(j+1)-((F1/RI)*(j-1)^2)*(-1)^(j);
       end
       rhs(bc_pos)=bc_vals(l)*F1^2;
   elseif strcmp(typ,'SF2_R')  %stress free for v at right v_r - (v)/r
       for j=1:M
           OP(bc_pos,j)=(j-1)^2-(F1/RO)*(1);  %RO
       end
       rhs(bc_pos)=bc_vals(l)*F1;
   elseif strcmp(typ,'SF2_L') %stress free for v at left
       for j=1:M
           OP(bc_pos,j)=((j-1)^2*(-1)^(j))-((F1/RI)*(-1)^(j+1));
       end
       rhs(bc_pos)=bc_vals(l)*F1;
   elseif strcmp(typ,'SF3_R')  %stress free for potential function at right phi_rrr + (phi_rr)/r - (phi_r)/r^2
       for j=1:M
           OP(bc_pos,j)=(1/15)*(j-1)^2*((j-1)^2-1)*((j-1)^2-4)+(F1/RO)*(1/3)*(j-1)^2*((j-1)^2-1)-(F1^2/RO^2)*((j-1)^2);  %RO
       end
       rhs(bc_pos)=bc_vals(l)*F1^3;
   elseif strcmp(typ,'SF3_L') %stress free for potential function at left
       for j=1:M
           OP(bc_pos,j)=(1/15)*((j-1)^2*((j-1)^2-1)*((j-1)^2-4))*(-1)^(j)+(F1/RI)*((1/3)*(j-1)^2*((j-1)^2-1))*(-1)^(j+1)-((F1^2/RI^2)*(j-1)^2)*(-1)^(j);
       end
       rhs(bc_pos)=bc_vals(l)*F1^3;
   elseif strcmp(typ,'SF4_R')
       for j=1:M
           OP(bc_pos,j)=(1/3)*(j-1)^2*((j-1)^2-1)+(F1/RO)*(j-1)^2-(F1^2/RO^2)*1;
       end
       rhs(bc_pos)=bc_vals(l)*F1^2;
   elseif strcmp(typ,'SF4_L')
       for j=1:M
           OP(bc_pos,j)=(1/3)*((j-1)^2*((j-1)^2-1))*(-1)^(j+1)+(F1/RI)*(j-1)^2*(-1)^(j)-(F1^2/RI^2)*(-1)^(j+1);
       end
       rhs(bc_pos)=bc_vals(l)*F1^2;
   elseif strcmp(typ,'D3_R') % third derivative
       for j=1:M
           OP(bc_pos,j)=(1/15)*(j-1)^2*((j-1)^2-1)*((j-1)^2-4);
       end
       rhs(bc_pos)=bc_vals(l)*F1^3;
   elseif strcmp(typ,'D3_L')
       for j=1:M
           OP(bc_pos,j)=(1/15)*((j-1)^2*((j-1)^2-1)*((j-1)^2-4))*(-1)^(j);
       end
       rhs(bc_pos)=bc_vals(l)*F1^3;
   elseif strcmp(typ,'I_RRR')
       for j=1:M
           OP(bc_pos,j)=RO*((1/15)*(j-1)^2*((j-1)^2-1)*((j-1)^2-4))-RI*((1/15)*((j-1)^2*((j-1)^2-1)*((j-1)^2-4))*(-1)^(j));
       end
       rhs(bc_pos)=bc_vals(l)*F1^3;
   elseif strcmp(typ,'intR') % Specify Integration in r in a cylinder
       %Cylindrical contribution
       for j=1:2:M-3
           OP(bc_pos,j+1)=(F1^2)*(-1)*((1/((j-2)*j))+(1/(j*(j+2))));
       end
       OP(bc_pos,M-1)=(F1^2)*(-1)*((1/((M-3)*(M-1)))+(1/(2*(M-1))));
       %Cartestian contribution
       for j=1:2:M-2
           OP(bc_pos,j)=2*F1*G1*(-1/(j*(j-2)));
       end
       OP(bc_pos,M)=2*F1*G1*(-1/(2*(M-1)));
       OP(bc_pos,1)=OP(bc_pos,1);
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'I_SF2_tilde')
       for j=1:M
           OP(bc_pos,j)=(1/RO)*(1)-(1/RI)*(-1)^(j+1);
       end
       OP(bc_pos,:)=OP(bc_pos,:)*1.0;
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'C1')
       % specify the value of C1
       OP(bc_pos,:)=zeros(1,M);
       OP(bc_pos,2)=1;
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'C0')
       % specify the value of C0
       OP(bc_pos,:)=zeros(1,M);
       OP(bc_pos,1)=1;
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'intV')  % specify the integral [r^2*f(r)]
       prod=zeros(1,M);
       % cartesian component
       for j=1:2:M
           prod(j)=2*(-1/(j*(j-2))); 
       end
       OP(bc_pos,:)=prod*(G1^2*F1);
       prod=zeros(1,M);
       % cylindrical component
       for j=1:2:M-1
          prod(j+1)=(-1)*((1/((j-2)*j))+(1/(j*(j+2)))); 
       end
       OP(bc_pos,:)=prod*(2*G1*F1^2)+OP(bc_pos,:);
       prod=zeros(1,M);
       % "r^2" component
       for j=0:2:M-1
          prod(j+1)=(-1)*((1/(2*(j-1)*(j-3)))+(1/((j-1)*(j+1)))+(1/(2*(j+1)*(j+3)))); 
       end
       OP(bc_pos,:)=prod*F1^3+OP(bc_pos,:);
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'ASYM')
       % Anti-Symmetric Top/Bottom
       for j=1:2:M
           OP(bc_pos,j)=2;
       end
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'SYM')
       % Anti-Symmetric Top/Bottom
       for j=2:2:M
           OP(bc_pos,j)=2;
       end
       rhs(bc_pos)=bc_vals(l);
   elseif strcmp(typ,'ZERO')
       % Leave a line intentially full of zeros
       % do nothing for this entry
   else
       disp(sprintf('Invalid bc_type at entry : %d',l))
   end
end

% OP(end-3:end,:)
% pause