function wg_dat = shape_fun(x,args)

% Shape function and coupling constant 
% Kuldeep Amarnath
% August, 2002
%
% Returns the shape function y = f(x), where x,y are the parametric coordinates
% the two waveguides [1]. Also calculates the local coupling coeff. using a fitted
% equation derived from LoCoup.m. This function is called by ToCoup.m.
%
% [1] Matsuhara, M., Watanabe, A., J. Opt. Soc. Amer., V65, n2, p 163, (1975)

%***** Racetrack shaped:  |C ******

% define the shape of waveguide using inputs or defaults

if length(args) == 7
    lambda = args{1};
    R = args{2};
    L = args{3};
    D = args{4};
    H = args{5};
    n1_eff = args{6};
    n2_eff = args{7};
else % use defaults
    lambda = 1.55e-06; % wavelength
    R = 10.0e-06; % radius of curved portion
    L = 0; % length of straight portion
    D = 0; % minimum horz. dist. between ring and bus.
    H = 0.75e-06; % vert. dist. between ring and bus.
    n1_eff = 3.08;
    n2_eff = 3.08;
end

% ****** end of user input *******

x0 = -2*R - D -L/2;
xN = 2*R + D +L/2;

% intialize the calling program
if ischar(x)
    beta_1 = 2*pi/lambda * n1_eff;
    beta_2 = 2*pi/lambda * n2_eff;
    wg_dat = [beta_1, beta_2,x0,xN,R,L];
    return;
end
        
% evaluate y,c1,c2
if x <= -L/2
    % eval. y = f(x)
    y = 2*R * atan((x+L/2)/(2*R + D)) - L/2;

    % distance between corresponding segments in x and y.
    x_coord = [(2*R+D)*tan((y+L/2)/(2*R)), 0];
    y_coord = [R*sin((y+L/2)/R), D + R*(1-cos((y+L/2)/R))];
    d = sqrt((x_coord(1)-y_coord(1))^2 + (x_coord(2)-y_coord(2))^2);
    
    % local coup. coeff. derived using LoCoup.m and OWMS waveguide mode simulator
    c0 = CoupFun(d,H,lambda);
    dydx = 2*R/(2*R+D) * (cos((y+L/2)/(2*R)))^2;
    c1 = c0 * cos((y+L/2)/R) * sqrt(dydx);
    c2 = c0 * cos((y+L/2)/R) * 1/sqrt(dydx);
    
elseif abs(x) < L/2
    y = x;
    
    d = D;
    
    c0 = CoupFun(d,H,lambda);
    dydx = 1;
    c1 = c0;
    c2 = c0;
    
else % x > L/2
    % y = f(x)
    y = 2*R * atan((x-L/2)/(2*R + D)) + L/2;     

    % distance between corresponding segments in x and y.
    x_coord = [(2*R+D)*tan((y-L/2)/(2*R)), 0];
    y_coord = [R*sin((y-L/2)/R), D + R*(1-cos((y-L/2)/R))];
    d = sqrt((x_coord(1)-y_coord(1))^2 + (x_coord(2)-y_coord(2))^2);
    
    % local coup. coeff. derived using LoCoup.m and OWMS waveguide mode simulator
    c0 = CoupFun(d,H,lambda);
    dydx = 2*R/(2*R+D) * (cos((y-L/2)/(2*R)))^2;
    c1 = c0 * cos((y-L/2)/R) * sqrt(dydx);
    c2 = c0 * cos((y-L/2)/R) * 1/sqrt(dydx);
       
end

wg_dat = [y,c1,c2];

% function to return the coupling coeff. C0 = f(x,y,lambda)
function CoupCoeff = CoupFun(x,y,lambda)

% lambda at which mode solver simulation was carried out
lambda0 = 1.55e-06;

% Output of LoCoup.m
CoupCoeff = 5.611e6 * exp(-5.316e6 * sqrt(x*x+y*y));