function CoupFun = LoCoup(file1,file2pos1,xoff1,yoff1,file2pos2,xoff2,yoff2)
% Calculation of Local Coupling Coefficients of Parallel Waveguides, CMD line
% Kuldeep Amarnath
% August, 2002
% 
% Takes in three OWMS optical mode solver output files (usually olink.dat) 
% and calculates the coupling coefficient as described by Okamoto [1]
% with a slightly different normalization of the electric field. The OWMS
% solutions for the two waveguides should be based on the same mesh. It
% then fits an exponential function, which is returned as a string.
%
% Ref: [1]  Katsunari Okamoto,'Fundamentals of Optical Waveguides',
%           Academic Press,2000. 

% define constants
eps0=8.84e-12;
mu0 = 4*pi*1e-07;
C = 2.997e08;

if nargin ~= 7
    error('Usage: LoCoup(''waveguide1.dat'',''waveguide2pos1.dat'',xoff1,yoff1,''waveguide2pos2.dat'',xoff2,yoff2)')
else
    fid1 = fopen(file1,'r');
    fid2 = fopen(file2pos1,'r');
    fid3 = fopen(file2pos2,'r');
end


% read header info. file1

lambda = fscanf(fid1,'%f // wavelength',1) * 1e-06;
line = fgets(fid1);
mesh = fscanf(fid1, '%f %f 1 // mesh');
nx = mesh(1);
ny = mesh(2);

for i=1:6
    line = fgets(fid1);
end

n_eff1 = fscanf(fid1,'%f',1);
line = fgets(fid1);

% read data 

for col=1:nx
    for row=1:ny
        data = fscanf(fid1,'%f',7);
        Xpos(col) = data(1);
        Ypos(row) = data(2);
        indx1(row,col) = data(4);
        E1_x(row,col) = data(6);
        E1_y(row,col) = data(7);
    end
end

fclose(fid1);

% do same for second file

% skip header
for i=1:7
    line = fgets(fid2);
end
n_eff2 = fscanf(fid2,'%f',1);
line = fgets(fid2);

% read data 
for col=1:nx
    for row=1:ny
        data = fscanf(fid2,'%f',7);
        indx2(row,col) = data(4);
        E2_x(row,col) = data(6);
        E2_y(row,col) = data(7);
    end
end
fclose(fid2);

% do same for third file

% skip header
for i=1:8
    line = fgets(fid3);
end

% read data 
for col=1:nx
    for row=1:ny
        data = fscanf(fid3,'%f',7);
        indx3(row,col) = data(4);
        E3_x(row,col) = data(6);
        E3_y(row,col) = data(7);
    end
end
fclose(fid3);

% calculate the coupling coefficient for wg_1 and wg_2_pos1

indx12 = max(indx1,indx2); % the index profile with both waveguides present

I = 0;
P1 = 0;
P2 = 0;
for row=2:ny-1
    for col=2:nx-1
        dx = abs(Xpos(col+1)-Xpos(col-1))/2;
        dy = abs(Ypos(row+1)-Ypos(row-1))/2;
        dI = (indx12(row,col)^2 - indx2(row,col)^2) * ...
             (E1_x(row,col)*E2_x(row,col) + E1_y(row,col)*E2_y(row,col)) * dx*dy;
        I = I + dI;
        % normalization constants (Power flow in Z direction)
        dP1 = (E1_x(row,col)*E1_x(row,col) + E1_y(row,col)*E1_y(row,col))*dx*dy;
        dP2 = (E2_x(row,col)*E2_x(row,col) + E2_y(row,col)*E2_y(row,col))*dx*dy;
        P1 = P1 + dP1;
        P2 = P2 + dP2;
    end
end

P1 = 1/2 * n_eff1/(C*mu0) * P1;
P2 = 1/2 * n_eff2/(C*mu0) * P2;
I = I/(sqrt(P1)*sqrt(P2));

CCoeff12 = 2.997e08*2*pi/lambda * eps0 * I/4;

% calculate the coupling coefficient for wg_1 and wg_2_pos2

indx13 = max(indx1,indx3); % the index profile with both waveguides present

I = 0;
P1 = 0;
P3 = 0;
for row=2:ny-1
    for col=2:nx-1
        dx = abs(Xpos(col+1)-Xpos(col-1))/2;
        dy = abs(Ypos(row+1)-Ypos(row-1))/2;
        dI = (indx13(row,col)^2 - indx3(row,col)^2) * ...
             (E1_x(row,col)*E3_x(row,col) + E1_y(row,col)*E3_y(row,col)) * dx*dy;
        I = I + dI;
        % normalization constants (Power flow in Z direction)
        dP1 = (E1_x(row,col)*E1_x(row,col) + E1_y(row,col)*E1_y(row,col))*dx*dy;
        dP3 = (E3_x(row,col)*E3_x(row,col) + E3_y(row,col)*E3_y(row,col))*dx*dy;
        P1 = P1 + dP1;
        P3 = P3 + dP3;
    end
end

P1 = 1/2 * n_eff1/(C*mu0) * P1;
P3 = 1/2 * n_eff2/(C*mu0) * P3;
I = I/(sqrt(P1)*sqrt(P3));

CCoeff13 = 2.997e08*2*pi/lambda * eps0 * I/4;

% fit coupling function C = C0 * exp(-C1*d)

d1 = sqrt(xoff1*xoff1+yoff1*yoff1);
d2 = sqrt(xoff2*xoff2+yoff2*yoff2);

C1 = 1/(d2-d1) * log(CCoeff12/CCoeff13);
C0 = CCoeff12 * exp(C1*d1);

CoupFun = [num2str(C0,'%0.3e'),' * exp(-',num2str(C1,'%0.3e'),' * sqrt(x*x + y*y))'];