## Generate a spiral line ##
function mx_poly_spiral(r,theta,coord,order,para){
    #UNTITLED2 Summary of this function goes here
    #   Detailed explanation goes here

    dL = para.dL;
    r_init = r(1);
    r_end = r(2);
    theta_init = theta(1);
    theta_end = theta(2);
    x_init = coord(1);
    y_init = coord(2);
    
    K0 = 1/r_init;
    K1 = 1/r_end;
    
    L0 = abs(theta_end - theta_init)/(K0+(K1-K0)*order/(order+1));
    
    L = [0:dL:L0];
    K = K0 + (K1 - K0)/L0^order * (L0^order - abs(L-L0)^order);
    
    R = 1/K;
    R = (R<=para.R_max)*R + (R>para.R_max)*para.R_max*ones(length(R),1);
    direction = sign(theta_end-theta_init);
    dt = direction*dL/R;
    #theta_temp = cumsum(dt)+theta_init;
    theta_temp = dt;
    
    x=zeros(length(L),1)+x_init;
    y=zeros(length(L),1)+y_init;
    
    
    
    for (i=2;i<=length(L);i=i+1){
        theta_temp(i) = theta_temp(i)+theta_temp(i-1);
        cur_theta = theta_temp(i)+theta_init;
        pre_theta = theta_temp(i-1)+theta_init;
        x(i) = x(i-1) + direction* R(i)*( sin( cur_theta ) - sin(pre_theta ) );
        y(i) = y(i-1) - direction* R(i)*( cos( cur_theta ) - cos( pre_theta ) );
    
    }
    
    theta_temp = [theta_temp(1);theta_temp(2:50:end-1);theta_temp(end)]+theta_init;
    x = [x(1);x(2:50:end-1);x(end)];
    y = [y(1);y(2:50:end-1);y(end)];
    
    
    vtx = [x,y,theta_temp];
    
    return vtx;    
    
}

function mx_wg_draw(vtx,width){ 
    #UNTITLED6 Summary of this function goes here
    #   Detailed explanation goes here
    
    z = vtx(:,1) + 1i*vtx(:,2); # complex points
    #dz = diff(z); # direction of each point
    dz = z(2:end) - z(1:end-1);    
    
    dz = [transpose(dz),dz(end)];
    dir_upper = -1i*real(dz)+imag(dz);
    dir_down =  1i*real(dz)-imag(dz);
    
    p_upper = [z + dir_upper*width/2/abs(dir_upper)];
    p_down = [z+ dir_down*width/2/abs(dir_down)];
    
    wg = struct;
    wg.curve_inner = [real(p_upper),imag(p_upper)];
    wg.curve_outer = [real(p_down),imag(p_down)]; 
    
    return wg;
    
}

function mx_euler_wg(vtx,width,offset){ 
    #UNTITLED6 Summary of this function goes here
    #   Detailed explanation goes here
    
    z = vtx(:,1) + 1i*vtx(:,2); # complex points
    #dz = diff(z); # direction of each point
    dz = sin(vtx(:,3))*1i + cos(vtx(:,3));   
    
    dz = [transpose(dz)];
    dir_upper = -1i*real(dz)+imag(dz);
    dir_down =  1i*real(dz)-imag(dz);
    
    p_upper = [z + dir_upper*(offset+width/2)/abs(dir_upper)];
    p_down = [z+ dir_down*(-offset+width/2)/abs(dir_down)];
    
    wg = struct;
    wg.curve_inner = [real(p_upper),imag(p_upper)];
    wg.curve_outer = [real(p_down),imag(p_down)]; 
    
    return wg;
    
}



function mx_euler_wg2wg(euler_para,bend_angle,theta_start,coord,bend_type,width_type,Height,Material,vtx_flip)
{
     R0 = euler_para.R0;
     R1 = euler_para.R1;
     Win = euler_para.w0;
     dW = euler_para.w1 - euler_para.w0;
     order = euler_para.order;
     
     if (bend_type=='single'){
        vtx_start = mx_poly_spiral([R0,R1],[theta_start,bend_angle+theta_start],[0,0],order,euler_para.para);
        p_start = vtx_start(1,:);
        p_end = vtx_start(end,:);
        vtx_euler_bend = vtx_start;
     }
        
     else {
        R2 = euler_para.R2;

        vtx_start = mx_poly_spiral([R0,R1],[theta_start,bend_angle/2+theta_start],[0,0],order,euler_para.para);
        vtx_stop = mx_poly_spiral([R1,R2],[bend_angle/2+theta_start,bend_angle+theta_start],[vtx_start(end,1),vtx_start(end,2)],order,euler_para.para);
        p_start = vtx_start(1,:);
        p_end = vtx_stop(end,:);
        vtx_euler_bend = [vtx_start;vtx_stop(2:end,:)] ;
    }
    
    ## attaching waveguide    
    dx = abs(p_end(2) - p_start(1));## displacement in x direction
    dL = (vtx_euler_bend(2:end,2) - vtx_euler_bend(1:end-1,2))^2 + (vtx_euler_bend(2:end,1) - vtx_euler_bend(1:end-1,1))^2;
    dL = sqrt(dL);
    ##L = cumsum(dL) ## L for each pieces
    L = zeros(length(dL),1);
    L(1) = dL(1);
    for (idx=2;idx<=length(L);idx=idx+1){
      L(idx) = L(idx-1)+dL(idx);            
        }        
    L = [0;L];

    L0 = sum(dL);
    w_offset = 0; 
        if (width_type=='cos'){## in this situation, dW is the difference of input and output
            dy = abs(p_end(2) - p_start(2)); ## displacement in y direction
            vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2);
            z = vtx_euler_bend(:,3);

            #z = vtx_euler_bend(:,1) + 1i*vtx_euler_bend(:,2);
            w = dW/2*cos(z*pi/abs(bend_angle)) + (Win*2+dW)/2;
        }
        
        else if (width_type=='sin'){ ## in this situation, win = wout, dW is the middle width difference
            if (abs(bend_angle-pi)<0.001){                
                dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction
            }
            else {
            dy = abs(p_end(2) - p_start(2)); ## displacement in y direction                    
                }                
            vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2);
            z = vtx_euler_bend(:,3);
            
            w = dW*cos(z+pi/2)^2 + Win; ## revised 2023.03.27
            
            vtx_euler_bend(:,2) =  dy + vtx_euler_bend(:,2);

            }

        else if (width_type=='pumpkin'){ ## in this situation, win = wout, dW is the middle width difference
            if (abs(bend_angle-pi)<0.001){                
                dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction
            }
            else {
            dy = abs(p_end(2) - p_start(2)); ## displacement in y direction                    
                }                
            vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2);
            z = vtx_euler_bend(:,3);
            
            z = z;    
            z = z^0.5*(pi/2)^0.5;     
            z = sin(z)^2*pi/2;     
            
            w = dW*sin( z )^2 + Win; ## revised 2023.05.04
            #w = dW*sin( z )^2 + Win; ## revised 2023.05.04
            #w = dW*theta/(pi/2) + Win; ## revised 2023.05.04
            vtx_euler_bend(:,2) =  dy + vtx_euler_bend(:,2);

            }

        else if (width_type=='special'){ ## in this situation, win = wout, dW is the middle width difference
            if (abs(bend_angle-pi)<0.001){                
                dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction
            }
            else {
            dy = abs(p_end(2) - p_start(2)); ## displacement in y direction                    
                }                
            vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2);
            z = vtx_euler_bend(:,3);
            
            z = z;    
            z = z^0.65*(pi/2)^0.35;     
            z = sin(z)^2*pi/2;     
            
            w = dW*sin( z )^2 + Win; ## revised 2023.05.04
            vtx_euler_bend(:,2) =  dy + vtx_euler_bend(:,2);

            }
            
        else if (width_type=='sin2'){ ## in this situation, win = wout, dW is the middle width difference
            if (abs(bend_angle-pi)<0.001){                
                dy = abs(vtx_start(end,2) - vtx_start(1,2)); ## displacement in y direction
            }
            else {
            dy = abs(p_end(2) - p_start(2)); ## displacement in y direction                    
                }                
            vtx_euler_bend(:,2) = -dy + vtx_euler_bend(:,2);
            z = vtx_euler_bend(:,1) + 1i*vtx_euler_bend(:,2);
            w = dW/2*sin(angle(z)*2+abs(bend_angle))^2 + Win+dW/2;
            }

        else if (width_type=='linear'){
            w = dW/L0*L + Win;}

        else if (width_type=='dual_linear'){
            w = dW/2/(L0/2)*abs(L-L0/2) + Win+dW/2;}
            
        else if (width_type=='linear_offset'){
            w = dW/L0*L + Win;
            w_offset = euler_para.w_offset; 
            }

        else { ## default linear from input to output
            w = dW/L0*L + Win;}

          
        sz = abs([p_end(1) - p_start(1),p_end(2) - p_start(2)]); ## the size of the bending
        wg = mx_euler_wg(vtx_euler_bend,w,w_offset);
        
        vtx = [wg.curve_outer;flip(wg.curve_inner,1)];
        
        if (vtx_flip(1)==1){
            vtx(:,1) = -vtx(:,1);}
        if (vtx_flip(2)==1){
            vtx(:,2) = -vtx(:,2);}
            
        mx_poly('euler',coord,vtx,Height,Material,1,0);  
        wg = struct;
        wg.sz = sz;
        wg.w = w;
        wg.vtx = vtx_euler_bend;  ## central line
        
        z = vtx(:,1) + 1i*vtx(:,2); # complex points
        #dz = diff(z); # direction of each point
        dz = z(2:end) - z(1:end-1);    
        
        dz = [transpose(dz),dz(end)];
        dir_upper = -1i*real(dz)+imag(dz);
        dir_down =  1i*real(dz)-imag(dz);
    
        wg.angle = -angle(dir_upper);
        return wg;
    
}