"""Curve and spiral geometry primitives."""

from typing import Any, Optional
from cmath import pi

import nazca as nd
import numpy as np

from .polygons import _my_polygon

class Conchoid:
    def __init__(self,R0: Any,kR: Any,T: Any,w: float,layer: str,w_end: Optional[float]=None,res: float=0.1,final_flat: Any=None,begin_flat: Any=None,xs: Optional[str]=None) -> None:
        ## with half circle to be one cycle
        
        if (w_end==None):
            w_end = w
            
        
        with nd.Cell(instantiate=False)as C:
            n_sects = int(np.floor(T/np.pi)) ## intersecting into different semi-circle 
            
            L = R0*T + 1/2*kR*T**2  ## The total length of the Conchoid center line 
            n_points = int(np.floor(L/res))+1

            ## calculating sections
            if (np.abs(T-n_sects*np.pi) <0.0001):
                n_sects = n_sects
            else:
                n_sects = n_sects+1

            # res_sect = int(np.floor(res/n_sects))+1
            
            # res = 0

            if (layer!=None):
                nd.add_xsection(name='temp')
                nd.add_layer2xsection(xsection='temp',layer=layer,growx=0,growy=0)
                xs = 'temp'

            for _n_ in range(0,n_sects):
                # phi_start = _n_*pi
                # phi_end = min(T,_n_*pi+pi) ## forward placement
                
                phi_end = T - _n_*pi
                phi_start = max(0,T - _n_*pi - pi)                
                
                L_sect = R0*(phi_end-phi_start) + 1/2*kR*(phi_end**2 - phi_start**2)
                n_points = int(L_sect/res)+1
                
                if (layer!=None):
                    Theta = np.linspace(phi_start, phi_end ,n_points)
                    R = (Theta*kR+R0) ## conchoid function 
                    
                    # res = kR/2*T*T + R0*T
                    # res = res + np.sum(R[0:-1]*np.diff(Theta))
                    
                    w_cur = w
                    
                    if (_n_==0):
                        w_cur = np.linspace(w,w_end,n_points)
                        # print("Loading Taper area")

                    vtx_cx = R*np.cos(Theta)
                    vtx_cy = R*np.sin(Theta)
                    vtx_center =  np.c_[vtx_cx,vtx_cy]
                    
                    e_theta = -1/((R0/kR)+Theta)  ## actuall norm towards spiral
                    # e_rou = np.ones(len(e_theta))
                    ey = np.sin(Theta) - np.cos(Theta)*kR/R
                    ex = np.cos(Theta) + np.sin(Theta)*kR/R
                    
                    if (final_flat!=None and _n_==0):
                        ey[-1] = np.sin(final_flat/180*np.pi) 
                        ex[-1] = np.cos(final_flat/180*np.pi) 
                        
                    if (begin_flat!=None and _n_==n_sects-1):
                        ey[0] = np.sin(begin_flat/180*np.pi) 
                        ex[0] = np.cos(begin_flat/180*np.pi)                   
                    # if (final_flat!=None and _n_==0):
                    #     e_theta[-1] = final_flat 
                    # if (begin_flat!=None and _n_==n_sects-1):
                    #     e_theta[0] = begin_flat 
                    
                    if (_n_==0):
                        self.Atilt = np.arctan(ey[0]/ex[0])/np.pi*180
                        # print("Atilt_conchoid = %.5f" % self.Atilt)
                    
                    Lnorm = np.sqrt(np.power(ex,2)+np.power(ey,2))
                    
                    vtx_x = R*np.cos(Theta)
                    vtx_y = R*np.sin(Theta)
                    
                    vtx_out_x = vtx_x + w_cur/2*ex/Lnorm
                    vtx_out_y = vtx_y + w_cur/2*ey/Lnorm

                    vtx_in_x = vtx_x - w_cur/2*ex/Lnorm
                    vtx_in_y = vtx_y - w_cur/2*ey/Lnorm            

                    vtx_in = np.c_[np.flip(vtx_in_x),np.flip(vtx_in_y)]
                    vtx_out = np.c_[vtx_out_x,vtx_out_y]
                    vtx = np.r_[vtx_out,vtx_in]   
                    _my_polygon(layer_wg=layer,vtx=vtx).put(0,0,0)       
                
                elif(layer==None and xs!=None):
                    for layers,growx,growy,acc in nd.layeriter(xs=xs):
                        (a1,b1), (a2,b2),c1,c2 = growx         
                        Theta = np.linspace(phi_start, phi_end ,n_points)
                        R = (Theta*kR+R0)
                        # res = kR/2*T*T + R0*T
                        # res = res + np.sum(R[0:-1]*np.diff(Theta)
                        
                        w_cur = w*(a1-a2) + (b1-b2)
                       
                        if (_n_==0):
                            w_cur = np.linspace(w,w_end,n_points)
                            w_cur = w_cur*(a1-a2) + (b1-b2)
                        
                        vtx_cx = R*np.cos(Theta)
                        vtx_cy = R*np.sin(Theta)
                        vtx_center =  np.c_[vtx_cx,vtx_cy]
                        
                        # e_theta = -1/((R0/kR)+Theta)
                        # e_rou = np.ones(len(e_theta))

                        # if (final_flat!=None and _n_==0):
                        #     e_theta[-1] = final_flat 
                        # if (begin_flat!=None and _n_==n_sects-1):
                        #     e_theta[0] = begin_flat 

                        
                        ey = np.sin(Theta)*R - np.sin(Theta)*kR
                        ex = np.cos(Theta)*R + np.cos(Theta)*kR
                        
                        if (final_flat!=None and _n_==0):
                            ey[-1] = np.sin(final_flat/180*np.pi) 
                            ex[-1] = np.cos(final_flat/180*np.pi) 
                            
                        if (begin_flat!=None and _n_==n_sects-1):
                            ey[0] = np.sin(begin_flat/180*np.pi) 
                            ex[0] = np.cos(begin_flat/180*np.pi)  
                                                   
                        Lnorm = np.sqrt(np.power(ex,2)+np.power(ey,2))
                        vtx_x = R*np.cos(Theta)
                        vtx_y = R*np.sin(Theta)
                        
                        vtx_out_x = vtx_x + w_cur/2*ex/Lnorm
                        vtx_out_y = vtx_y + w_cur/2*ey/Lnorm

                        vtx_in_x = vtx_x - w_cur/2*ex/Lnorm
                        vtx_in_y = vtx_y - w_cur/2*ey/Lnorm            

                        vtx_in = np.c_[np.flip(vtx_in_x),np.flip(vtx_in_y)]
                        vtx_out = np.c_[vtx_out_x,vtx_out_y]
                        vtx = np.r_[vtx_out,vtx_in]   
                        _my_polygon(layer_wg=layers,vtx=vtx).put(0,0,0)    
            
            Rmax = T*kR+R0          
            ## revised in 2026.06.07 by Qin Yue
            # legacy: nd.Pin(name="a1").put(R0,0,-90)
            nd.Pin(name="opt_a1",type="optical:").put(R0,0,-90)
            ## revised in 2026.06.07 by Qin Yue
            # legacy: nd.Pin(name="b1").put(Rmax*np.cos(T),Rmax*np.sin(T),(T/np.pi*180+90))
            nd.Pin(name="opt_b1",type="optical:").put(Rmax*np.cos(T),Rmax*np.sin(T),(T/np.pi*180+90))
            
        self.L = L
        self.cell =C
        self.vtx_center = vtx_center
        self.vtx = vtx
        self.K_end = (np.power(np.max(R),2) + 2*np.power(kR,2)) / np.power((np.power(np.max(R),2) + np.power(kR,2)),1.5)
        self.R_end = 1/self.K_end

def _line2wg_(x,y,wu,wd,theta,n_points):
    """ building waveguide with center line and side expansion

    Args:
        vtx_line (list[float]): the location of points, [x,y]
        width (list[float]): the expansion width of points vertical to the pointing vector, [wu,wd]
        theta (list[float]): the pointing angle, [theta], 0 represent right, 180 represent left
    """

    theta = theta*np.pi/180
    x_u = x+wu*np.cos(theta+pi/2)
    x_d = x+wd*np.cos(theta-pi/2)

    y_u = y+wu*np.sin(theta+pi/2)
    y_d = y+wd*np.sin(theta-pi/2)

    ### polygon section, reducing resolution
    sect = np.linspace(start= 0,stop= len(x_u)-1,num= n_points)
    sect = np.asarray(sect, dtype = int)

    x_u = x_u[sect]
    x_d = x_d[sect]
    y_u = y_u[sect]
    y_d = y_d[sect]

    vtx_u = np.c_[x_u,y_u]
    vtx_d = np.c_[x_d,y_d]

    

    vtx = np.r_[vtx_u,np.flip(vtx_d,0)]

    return vtx

def _my_poly_spiral(r,theta,order,res,R_max,sz_restrict=None):
        ''' generating a poly spiral curve
        Args
            r (2*1 list) :r[0] is the begining 
            theta (2*1 list) :theta[0] is the begining [in degree]
        Return
            frame (nazca.cell):
        '''
        theta[0] = theta[0]/180*np.pi ## angle format changing
        theta[1] = theta[1]/180*np.pi ## angle format changing

        K_ends = np.array([1/r[0],1/r[1]]) ## definition of the curvature, r[0] is the beginnin and r[1] is the ending
        L0 = np.abs(theta[0]-theta[1])/(K_ends[0] + (K_ends[1]-K_ends[0])*order/(order+1))
        L = np.linspace(0,L0,int(np.floor(L0/res)+1)) ## L = [0:res:L0];
        K = K_ends[0] + (K_ends[1] - K_ends[0])/np.power(L0,order)*(np.power(L0,order) - np.power(np.abs(L-L0),order))
        R = 1/K
        
        dir = np.sign(theta[1] - theta[0])
        dt = dir*res/R

        theta_temp = np.cumsum(dt) + theta[0]
        """ 2023.08.01 updated, using array calculation instead of for loop"""
        
        dx = dir*R[1:]*( np.sin(theta_temp[1:]) - np.sin(theta_temp[0:-1]))
        dy = -dir*R[1:]*( np.cos(theta_temp[1:]) - np.cos(theta_temp[0:-1]))
        
        x = np.r_[0,np.cumsum(dx)]
        y = np.r_[0,np.cumsum(dy)]
        # x = np.zeros(len(L))
        # y = np.zeros(len(L))
        
        # idx = np.linspace(1,len(L)-1,len(L)-1)

        # for _idx_ in idx :
        #     _idx_ = int(_idx_)
        #     x[_idx_] = x[_idx_-1] + dir*R[_idx_]*( np.sin(theta_temp[_idx_]) - np.sin(theta_temp[_idx_-1]))
        #     y[_idx_] = y[_idx_-1] - dir*R[_idx_]*( np.cos(theta_temp[_idx_]) - np.cos(theta_temp[_idx_-1]))

        vector = np.c_[x,y,theta_temp,L]

        return (vector,L0)

class Clothoid:
    def __init__(self,
                 name:str=None,
            R: 'list|np.ndarray'=[10,20],
            w: 'list|np.ndarray|float'=[0.4,0.5], ## w either has the length as R, or 1 element, or 2 element
            A: 'list|np.ndarray'=[0,45],
            width_type: str='sine',
            spiral_order: float=1,
            Rmax:float=10000,
            dL_cal:float=0.001,
            dL_wg: float = 0.1,
            # n_points:int=1024,
            xs:str='strip',
            layer:str=None,
            sharp_patch:bool=True,
            end_patch : bool=True,
            show_pins:bool=False) -> None:
        """_summary_

        Args:
            R (list|np.ndarray, optional): Curvature radius in each attaching point. Defaults to [10,20].
            w (list|np.ndarray|float, optional): Width at each attaching point corresponding to R, or it can be set to one or two element. Defaults to [0.4,0.5].
            A (list|np.ndarray, optional): Angle at each attaching point. Defaults to [0,45].
            width_type (str, optional): The width function with length or angle 'linear' 'linear2' 'sine' 'sine2'. Defaults to 'sine'.
            spiral_order (float, optional): The curvature order of spiral. Defaults to 1.
            Rmax (float, optional): Maxmum radius. Defaults to 10000.
            res (float, optional): Resolution in calculation. Defaults to 0.001.
            n_points (int, optional): Resolution in GDS. Defaults to 1024.
            xs (str, optional): XSection of the devices. Defaults to 'strip'.
            layer (str, optional): Layer of the devices. Defaults to None.
            sharp_patch (bool, optional): Either to patch. Defaults to True.
            show_pins (bool, optional): Either to show pins. Defaults to False.

        Raises:
            Exception: _description_
            Exception: _description_
        """

        if (isinstance(w,int) or isinstance(w,float)):
            w= np.array([w,w])   
            
        self.name = name     
        self.R = R
        self.A = A
        self.width_type = width_type
        self.spiral_order = spiral_order
        self.dL_cal = dL_cal
        # self.n_points = n_points
        self.xs =xs 
        self.layer = layer
        self.dL_wg = dL_wg
        
        if (len(R) != len(A)):
            raise Exception("ERROR: In <mxpic::structures::Colthoid>, <A> and <R> are not matched in length, please keep len(A) = len(R)")
        
        if (isinstance(spiral_order,int) or isinstance(spiral_order,float)):
            spiral_order = spiral_order*np.ones(len(R)-1)
            
        elif(isinstance(spiral_order,list)):
            spiral_order = np.array(spiral_order)
            
        ## center curve routing
        _idx_act_=0
        for _idx_ in range(0,len(R)-1):
            if ( abs( A[_idx_] - A[_idx_+1] )<0.001 ):

                continue
            vec_cur,L0_cur = _my_poly_spiral([R[_idx_],R[_idx_+1]],[A[_idx_],A[_idx_+1]],spiral_order[_idx_],dL_cal,Rmax) 
            _idx_act_ = _idx_act_+1

            x_cur = vec_cur[:,0]
            y_cur = vec_cur[:,1]
            theta_cur = vec_cur[:,2]/np.pi*180  ## pointing vector
            L_cur = vec_cur[:,3] 

            if (_idx_act_==1):
                L = L_cur
                x = x_cur
                y = y_cur
                theta = theta_cur
                L0 = L0_cur
            else :
                L = np.r_[L,L_cur+L[-1]]
                x = np.r_[x,x_cur+x[-1]]
                y = np.r_[y,y_cur+y[-1]]
                theta = np.r_[theta,theta_cur]     
                L0 = L0 + L0_cur           
            
            if (len(w)>2 and len(w)==len(R)):
                w_cur = (w[_idx_+1]-w[_idx_])/L0_cur*L_cur + (w[_idx_])
                if (_idx_act_==1):
                    w_fianl = w_cur
                else :
                    w_fianl = np.r_[w_fianl,w_cur]

        self.x = x
        self.y = y
        self.L = L
        self.L0 = L0
        self.theta = theta
        self.vtx_center = np.c_[x,y]   
        self.end_patch = end_patch
        self.sz = [np.abs(max(self.x) - min(self.x)),np.abs(max(self.y) - min(self.y))]
        
        if (dL_wg!=None):
            self.n_points = int(np.floor(self.L0/self.dL_wg)+1) ## overwrite n_points
        
        # print("n points",self.n_points)
        
        if (len(w)==2):
            ## width winding
            if (width_type=='linear'):
                w = (w[1]-w[0])/L0*L + w[0]
            elif (width_type=='dual_linear'):
                w = (w[1]-w[0])/L0/2*np.abs(L-L0/2) + w[0]
            elif (width_type=='sine'):
                w = (w[0]-w[1])*np.cos(theta/180*pi)*np.cos(theta/180*pi) + w[1]  
            elif (width_type=='dual_sine'):
                w = (w[0]-w[1])*np.cos(theta/2/180*pi)*np.cos(theta/2/180*pi) + w[1] 
            elif (width_type=='crow_customize' or width_type=='pumpkin'):
                dw = (w[1]-w[0])
                z = theta/180*np.pi
                z = np.sqrt(z)*np.sqrt(np.pi/2)
                z = np.sin(z)**2*np.pi/2

                w = dw*np.sin(z)**2 + w[0]

                
            else :
                w = (w[1]-w[0])/L0*L + w[0]     
            
            self.w = np.array(w) 
        elif (len(w)==len(R)):
            self.w = w_fianl
        else:
            raise Exception("ERROR, In <mxpic::structures::Clothoid>, <w> is not matched with <R>, please keep len(w)=2 or len(w)=R or w=int")
             
        self.cell = self.generate_gds(sharp_patch=sharp_patch,show_pins=show_pins)
  
        
    def generate_gds(self,sharp_patch,show_pins):
        if (self.name is None):
            self.instantiate = False
        else:
            self.instantiate = True
        with nd.Cell(name=self.name,instantiate=self.instantiate) as C:
            
            if (self.layer==None and self.xs!=None): ## if definition is in layers
                for layers,growx,growy,acc in nd.layeriter(xs=self.xs):
                    (a1,b1), (a2,b2),c1,c2 = growx
                    if (b1!=0 and b2!=0):
                        vtx_wg = _line2wg_(x=self.x,y=self.y,wu=self.w*a1+b1,wd= -self.w*a2-b2,theta=self.theta,n_points=self.n_points)

                        dX = np.max(vtx_wg[:,0]) - np.min(vtx_wg[:,0])
                        dY = np.max(vtx_wg[:,1]) - np.min(vtx_wg[:,1])
                        cX = np.max(vtx_wg[:,0])/2 + np.min(vtx_wg[:,0])/2
                        cY = np.max(vtx_wg[:,1])/2 + np.min(vtx_wg[:,1])/2
                        
                        if (sharp_patch):
                            if (self.end_patch):
                                nd.strt(length = dX+(b1-b2),width = dY,layer=layers).put(cX-dX/2-b1,cY,0)

                            else:
                                nd.strt(length = dX,width = dY,layer=layers).put(cX-dX/2,cY,0)     
                        else :
                            

                            _my_polygon(layers,vtx_wg).put(0,0,0)                   
                    
                    else : 
                        vtx_wg = _line2wg_(x=self.x,y=self.y,wu=self.w*a1+b1,wd= -self.w*a2-b2,theta=self.theta,n_points=self.n_points)
                        self.vtx =vtx_wg
                        _my_polygon(layers,vtx_wg).put(0,0,0)                   

                nd.Pin(name='a0',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
                nd.Pin(name='b0',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
                
                ## revised in 2026.06.07 by Qin Yue
                # legacy: nd.Pin(name='a1',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
                nd.Pin(name='opt_a1',width=self.w[0],type="optical:").put(self.x[0],self.y[0],self.A[0]+180)
                ## revised in 2026.06.07 by Qin Yue
                # legacy: nd.Pin(name='b1',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
                nd.Pin(name='opt_b1',width=self.w[-1],type="optical:").put(self.x[-1],self.y[-1],self.A[-1])
            elif(self.layer!=None) : ## if definition is in xsections

                vtx_wg = _line2wg_(x=self.x,y=self.y,wu=self.w/2,wd= self.w/2,theta=self.theta,n_points=self.n_points)

                _my_polygon(self.layer,vtx_wg).put(0,0,0)
                
                nd.Pin(name='a0',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
                nd.Pin(name='b0',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
                
                ## revised in 2026.06.07 by Qin Yue
                # legacy: nd.Pin(name='a1',width=self.w[0],type="Optical:").put(self.x[0],self.y[0],self.A[0]+180)
                nd.Pin(name='opt_a1',width=self.w[0],type="optical:").put(self.x[0],self.y[0],self.A[0]+180)
                ## revised in 2026.06.07 by Qin Yue
                # legacy: nd.Pin(name='b1',width=self.w[-1],type="Optical:").put(self.x[-1],self.y[-1],self.A[-1])
                nd.Pin(name='opt_b1',width=self.w[-1],type="optical:").put(self.x[-1],self.y[-1],self.A[-1])
                
            else:
                raise Exception("ERROR: In <mxpic::structures::Colthoid>, <layer | xs> not defined")
            
            if (show_pins):
                nd.put_stub()
            self.sz_p2p = [np.abs(self.x[-1] - self.x[0]),np.abs(self.y[-1] - self.y[0])]
            
        return C