"""
This module can be used to solve probelsm related to 2D parabolic arches
"""
from sympy.core.sympify import sympify
from sympy.core.symbol import Symbol, symbols
from sympy import diff, sqrt, cos, sin, atan, rad, Min
from sympy.core.relational import Eq
from sympy.solvers.solvers import solve
from sympy.functions import Piecewise
from sympy.plotting import plot
from sympy import limit
from sympy.utilities.decorator import doctest_depends_on
from sympy.external.importtools import import_module

numpy = import_module('numpy', import_kwargs={'fromlist':['arange']})

class Arch:
    """
    This class is used to solve problems related to a three hinged arch(determinate) structure.\n
    An arch is a curved vertical structure spanning an open space underneath it.\n
    Arches can be used to reduce the bending moments in long-span structures.\n

    Arches are used in structural engineering(over windows, door and even bridges)\n
    because they can support a very large mass placed on top of them.

    Example
    ========
    >>> from sympy.physics.continuum_mechanics.arch import Arch
    >>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
    >>> a.get_shape_eqn
    5 - (x - 5)**2/5

    >>> from sympy.physics.continuum_mechanics.arch import Arch
    >>> a = Arch((0,0),(10,1),crown_x=6)
    >>> a.get_shape_eqn
    9/5 - (x - 6)**2/20
    """
    def __init__(self,left_support,right_support,**kwargs):
        self._shape_eqn = None
        self._left_support  = (sympify(left_support[0]),sympify(left_support[1]))
        self._right_support  = (sympify(right_support[0]),sympify(right_support[1]))
        self._crown_x = None
        self._crown_y = None
        if 'crown_x' in kwargs:
            self._crown_x = sympify(kwargs['crown_x'])
        if 'crown_y' in kwargs:
            self._crown_y = sympify(kwargs['crown_y'])
        self._shape_eqn = self.get_shape_eqn
        self._conc_loads = {}
        self._distributed_loads = {}
        self._loads = {'concentrated': self._conc_loads, 'distributed':self._distributed_loads}
        self._loads_applied = {}
        self._supports = {'left':'hinge', 'right':'hinge'}
        self._member = None
        self._member_force = None
        self._reaction_force = {Symbol('R_A_x'):0, Symbol('R_A_y'):0, Symbol('R_B_x'):0, Symbol('R_B_y'):0}
        self._points_disc_x = set()
        self._points_disc_y = set()
        self._moment_x  = {}
        self._moment_y = {}
        self._load_x = {}
        self._load_y = {}
        self._moment_x_func = Piecewise((0,True))
        self._moment_y_func = Piecewise((0,True))
        self._load_x_func = Piecewise((0,True))
        self._load_y_func = Piecewise((0,True))
        self._bending_moment = None
        self._shear_force = None
        self._axial_force = None
        # self._crown = (sympify(crown[0]),sympify(crown[1]))

    @property
    def get_shape_eqn(self):
        "returns the equation of the shape of arch developed"
        if self._shape_eqn:
            return self._shape_eqn

        x,y,c = symbols('x y c')
        a = Symbol('a',positive=False)
        if self._crown_x and self._crown_y:
            x0  = self._crown_x
            y0 = self._crown_y
            parabola_eqn = a*(x-x0)**2 + y0 - y
            eq1 = parabola_eqn.subs({x:self._left_support[0], y:self._left_support[1]})
            solution = solve((eq1),(a))
            parabola_eqn = solution[0]*(x-x0)**2 + y0
            if(parabola_eqn.subs({x:self._right_support[0]}) != self._right_support[1]):
                raise ValueError("provided coordinates of crown and supports are not consistent with parabolic arch")

        elif self._crown_x:
            x0  = self._crown_x
            parabola_eqn = a*(x-x0)**2 + c - y
            eq1 = parabola_eqn.subs({x:self._left_support[0], y:self._left_support[1]})
            eq2 = parabola_eqn.subs({x:self._right_support[0], y:self._right_support[1]})
            solution = solve((eq1,eq2),(a,c))
            if len(solution) <2 or solution[a] == 0:
                raise ValueError("parabolic arch cannot be constructed with the provided coordinates, try providing crown_y")
            parabola_eqn = solution[a]*(x-x0)**2+ solution[c]
            self._crown_y = solution[c]

        else:
            raise KeyError("please provide crown_x to construct arch")

        return parabola_eqn

    @property
    def get_loads(self):
        """
        return the position of the applied load and angle (for concentrated loads)
        """
        return self._loads

    @property
    def supports(self):
        """
        Returns the type of support
        """
        return self._supports

    @property
    def left_support(self):
        """
        Returns the position of the left support.
        """
        return self._left_support

    @property
    def right_support(self):
        """
        Returns the position of the right support.
        """
        return self._right_support

    @property
    def reaction_force(self):
        """
        return the reaction forces generated
        """
        return self._reaction_force

    def apply_load(self,order,label,start,mag,end=None,angle=None):
        """
        This method adds load to the Arch.

        Parameters
        ==========

            order : Integer
                Order of the applied load.

                    - For point/concentrated loads, order = -1
                    - For distributed load, order = 0

            label : String or Symbol
                The label of the load
                - should not use 'A' or 'B' as it is used for supports.

            start : Float

                    - For concentrated/point loads, start is the x coordinate
                    - For distributed loads, start is the starting position of distributed load

            mag : Sympifyable
                Magnitude of the applied load. Must be positive

            end : Float
                Required for distributed loads

                    - For concentrated/point load , end is None(may not be given)
                    - For distributed loads, end is the end position of distributed load

            angle: Sympifyable
                The angle in degrees, the load vector makes with the horizontal
                in the counter-clockwise direction.

        Examples
        ========
        For applying distributed load

        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
        >>> a.apply_load(0,'C',start=3,end=5,mag=-10)

        For applying point/concentrated_loads

        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
        >>> a.apply_load(-1,'C',start=2,mag=15,angle=45)

        """
        y = Symbol('y')
        x = Symbol('x')
        x0 = Symbol('x0')
        # y0 = Symbol('y0')
        order= sympify(order)
        mag = sympify(mag)
        angle = sympify(angle)

        if label in self._loads_applied:
            raise ValueError("load with the given label already exists")

        if label in ['A','B']:
            raise ValueError("cannot use the given label, reserved for supports")

        if order == 0:
            if end is None or end<start:
                raise KeyError("provide end greater than start")

            self._distributed_loads[label] = {'start':start, 'end':end, 'f_y': mag}
            self._points_disc_y.add(start)

            if start in self._moment_y:
                self._moment_y[start] -= mag*(Min(x,end)-start)*(x0-(start+(Min(x,end)))/2)
                self._load_y[start] += mag*(Min(end,x)-start)
            else:
                self._moment_y[start] = -mag*(Min(x,end)-start)*(x0-(start+(Min(x,end)))/2)
                self._load_y[start] = mag*(Min(end,x)-start)

            self._loads_applied[label] = 'distributed'

        if order == -1:

            if angle is None:
                raise TypeError("please provide direction of force")
            height = self._shape_eqn.subs({'x':start})

            self._conc_loads[label] = {'x':start, 'y':height, 'f_x':mag*cos(rad(angle)), 'f_y': mag*sin(rad(angle)), 'mag':mag, 'angle':angle}
            self._points_disc_x.add(start)
            self._points_disc_y.add(start)

            if start in self._moment_x:
                self._moment_x[start] += self._conc_loads[label]['f_x']*(y-self._conc_loads[label]['y'])
                self._load_x[start] += self._conc_loads[label]['f_x']
            else:
                self._moment_x[start] = self._conc_loads[label]['f_x']*(y-self._conc_loads[label]['y'])
                self._load_x[start] = self._conc_loads[label]['f_x']

            if start in self._moment_y:
                self._moment_y[start] -= self._conc_loads[label]['f_y']*(x0-start)
                self._load_y[start] += self._conc_loads[label]['f_y']
            else:
                self._moment_y[start] = -self._conc_loads[label]['f_y']*(x0-start)
                self._load_y[start] = self._conc_loads[label]['f_y']

            self._loads_applied[label] = 'concentrated'


    def remove_load(self,label):
        """
        This methods removes the load applied to the arch

        Parameters
        ==========

        label : String or Symbol
            The label of the applied load

        Examples
        ========

        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
        >>> a.apply_load(0,'C',start=3,end=5,mag=-10)
        >>> a.remove_load('C')
        removed load C: {'start': 3, 'end': 5, 'f_y': -10}
        """
        y = Symbol('y')
        x = Symbol('x')
        x0 = Symbol('x0')

        if label in self._distributed_loads :

            self._loads_applied.pop(label)
            start = self._distributed_loads[label]['start']
            end = self._distributed_loads[label]['end']
            mag  = self._distributed_loads[label]['f_y']
            self._points_disc_y.remove(start)
            self._load_y[start] -= mag*(Min(x,end)-start)
            self._moment_y[start] += mag*(Min(x,end)-start)*(x0-(start+(Min(x,end)))/2)
            val = self._distributed_loads.pop(label)
            print(f"removed load {label}: {val}")

        elif label in self._conc_loads :

            self._loads_applied.pop(label)
            start = self._conc_loads[label]['x']
            self._points_disc_x.remove(start)
            self._points_disc_y.remove(start)
            self._moment_y[start] += self._conc_loads[label]['f_y']*(x0-start)
            self._moment_x[start] -= self._conc_loads[label]['f_x']*(y-self._conc_loads[label]['y'])
            self._load_x[start] -= self._conc_loads[label]['f_x']
            self._load_y[start] -= self._conc_loads[label]['f_y']
            val = self._conc_loads.pop(label)
            print(f"removed load {label}: {val}")

        else :
            raise ValueError("label not found")

    def change_support_position(self, left_support=None, right_support=None):
        """
        Change position of supports.
        If not provided , defaults to the old value.
        Parameters
        ==========

            left_support: tuple (x, y)
                x: float
                    x-coordinate value of the left_support

                y: float
                    y-coordinate value of the left_support

            right_support: tuple (x, y)
                x: float
                    x-coordinate value of the right_support

                y: float
                    y-coordinate value of the right_support
        """
        if left_support is not None:
            self._left_support = (left_support[0],left_support[1])

        if right_support is not None:
            self._right_support = (right_support[0],right_support[1])

        self._shape_eqn = None
        self._shape_eqn = self.get_shape_eqn

    def change_crown_position(self,crown_x=None,crown_y=None):
        """
        Change the position of the crown/hinge of the arch

        Parameters
        ==========

            crown_x: Float
                The x coordinate of the position of the hinge
                - if not provided, defaults to old value

            crown_y: Float
                The y coordinate of the position of the hinge
                - if not provided defaults to None
        """
        self._crown_x = crown_x
        self._crown_y = crown_y
        self._shape_eqn = None
        self._shape_eqn = self.get_shape_eqn

    def change_support_type(self,left_support=None,right_support=None):
        """
        Add the type for support at each end.
        Can use roller or hinge support at each end.

        Parameters
        ==========

            left_support, right_support : string
                Type of support at respective end

                    - For roller support , left_support/right_support = "roller"
                    - For hinged support, left_support/right_support = "hinge"
                    - defaults to hinge if value not provided

        Examples
        ========

        For applying roller support at right end

        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
        >>> a.change_support_type(right_support="roller")

        """
        support_types = ['roller','hinge']
        if left_support:
            if left_support not in support_types:
                raise ValueError("supports must only be roller or hinge")

            self._supports['left'] = left_support

        if right_support:
            if right_support not in support_types:
                raise ValueError("supports must only be roller or hinge")

            self._supports['right'] = right_support

    def add_member(self,y):
        """
        This method adds a member/rod at a particular height y.
        A rod is used for stability of the structure in case of a roller support.
        """
        if y>self._crown_y or y<min(self._left_support[1],  self._right_support[1]):
            raise ValueError(f"position of support must be between y={min(self._left_support[1],  self._right_support[1])} and y={self._crown_y}")
        x = Symbol('x')
        a = diff(self._shape_eqn,x).subs(x,self._crown_x+1)/2
        x_diff = sqrt((y - self._crown_y)/a)
        x1 = self._crown_x + x_diff
        x2 = self._crown_x - x_diff
        self._member = (x1,x2,y)

    def shear_force_at(self, pos = None, **kwargs):
        """
        return the shear at some x-coordinates
        if no x value provided, returns the formula
        """
        if pos is None:
            return self._shear_force
        else:
            x = Symbol('x')
            if 'dir' in kwargs:
                dir = kwargs['dir']
                return limit(self._shear_force,x,pos,dir=dir)
            return self._shear_force.subs(x,pos)

    def bending_moment_at(self, pos = None, **kwargs):
        """
        return the bending moment at some x-coordinates
        if no x value provided, returns the formula
        """
        if pos is None:
            return self._bending_moment
        else:
            x0 = Symbol('x0')
            if 'dir' in kwargs:
                dir = kwargs['dir']
                return limit(self._bending_moment,x0,pos,dir=dir)
            return self._bending_moment.subs(x0,pos)


    def axial_force_at(self,pos = None, **kwargs):
        """
        return the axial/normal force generated at some x-coordinate
        if no x value provided, returns the formula
        """
        if pos is None:
            return self._axial_force
        else:
            x = Symbol('x')
            if 'dir' in kwargs:
                dir = kwargs['dir']
                return limit(self._axial_force,x,pos,dir=dir)
            return self._axial_force.subs(x,pos)

    def solve(self):
        """
        This method solves for the reaction forces generated at the supports,\n
        and bending moment and generated in the arch and tension produced in the member if used.

        Examples
        ========

        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(10,0),crown_x=5,crown_y=5)
        >>> a.apply_load(0,'C',start=3,end=5,mag=-10)
        >>> a.solve()
        >>> a.reaction_force
        {R_A_x: 8, R_A_y: 12, R_B_x: -8, R_B_y: 8}

        >>> from sympy import Symbol
        >>> t = Symbol('t')
        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(16,0),crown_x=8,crown_y=5)
        >>> a.apply_load(0,'C',start=3,end=5,mag=t)
        >>> a.solve()
        >>> a.reaction_force
        {R_A_x: -4*t/5, R_A_y: -3*t/2, R_B_x: 4*t/5, R_B_y: -t/2}

        >>> a.bending_moment_at(4)
        -5*t/2
        """
        y = Symbol('y')
        x = Symbol('x')
        x0 = Symbol('x0')

        discontinuity_points_x = sorted(self._points_disc_x)
        discontinuity_points_y = sorted(self._points_disc_y)

        self._moment_x_func = Piecewise((0,True))
        self._moment_y_func = Piecewise((0,True))

        self._load_x_func = Piecewise((0,True))
        self._load_y_func = Piecewise((0,True))

        accumulated_x_moment = 0
        accumulated_y_moment = 0

        accumulated_x_load = 0
        accumulated_y_load = 0

        for point in discontinuity_points_x:
            cond = (x >= point)
            accumulated_x_load += self._load_x[point]
            accumulated_x_moment += self._moment_x[point]
            self._load_x_func = Piecewise((accumulated_x_load,cond),(self._load_x_func,True))
            self._moment_x_func = Piecewise((accumulated_x_moment,cond),(self._moment_x_func,True))

        for point in discontinuity_points_y:
            cond = (x >= point)
            accumulated_y_moment += self._moment_y[point]
            accumulated_y_load += self._load_y[point]
            self._load_y_func = Piecewise((accumulated_y_load,cond),(self._load_y_func,True))
            self._moment_y_func = Piecewise((accumulated_y_moment,cond),(self._moment_y_func,True))

        moment_A = self._moment_y_func.subs(x,self._right_support[0]).subs(x0,self._left_support[0]) +\
                   self._moment_x_func.subs(x,self._right_support[0]).subs(y,self._left_support[1])

        moment_hinge_left = self._moment_y_func.subs(x,self._crown_x).subs(x0,self._crown_x) +\
                            self._moment_x_func.subs(x,self._crown_x).subs(y,self._crown_y)

        moment_hinge_right = self._moment_y_func.subs(x,self._right_support[0]).subs(x0,self._crown_x)- \
                             self._moment_y_func.subs(x,self._crown_x).subs(x0,self._crown_x) +\
                             self._moment_x_func.subs(x,self._right_support[0]).subs(y,self._crown_y) -\
                             self._moment_x_func.subs(x,self._crown_x).subs(y,self._crown_y)

        net_x = self._load_x_func.subs(x,self._right_support[0])
        net_y = self._load_y_func.subs(x,self._right_support[0])

        if (self._supports['left']=='roller' or self._supports['right']=='roller') and not self._member:
            print("member must be added if any of the supports is roller")
            return

        R_A_x, R_A_y, R_B_x, R_B_y, T = symbols('R_A_x R_A_y R_B_x R_B_y T')

        if self._supports['left'] == 'roller' and self._supports['right'] == 'roller':

            if self._member[2]>=max(self._left_support[1],self._right_support[1]):

                if net_x!=0:
                    raise ValueError("net force in x direction not possible under the specified conditions")

                else:
                    eq1 = Eq(R_A_x ,0)
                    eq2 = Eq(R_B_x, 0)
                    eq3 = Eq(R_A_y + R_B_y + net_y,0)

                    eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
                             R_B_x*(self._right_support[1]-self._left_support[1])+moment_A,0)

                    eq5 = Eq(moment_hinge_right + R_B_y*(self._right_support[0]-self._crown_x) +\
                             T*(self._member[2]-self._crown_y),0)
                    solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

            elif self._member[2]>=self._left_support[1]:
                eq1 = Eq(R_A_x ,0)
                eq2 = Eq(R_B_x, 0)
                eq3 = Eq(R_A_y + R_B_y + net_y,0)
                eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
                         T*(self._member[2]-self._left_support[1])+moment_A,0)
                eq5 = Eq(T+net_x,0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

            elif self._member[2]>=self._right_support[1]:
                eq1 = Eq(R_A_x ,0)
                eq2 = Eq(R_B_x, 0)
                eq3 = Eq(R_A_y + R_B_y + net_y,0)
                eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])+\
                         T*(self._member[2]-self._left_support[1])+moment_A,0)
                eq5 = Eq(T-net_x,0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

        elif self._supports['left'] == 'roller':
            if self._member[2]>=max(self._left_support[1], self._right_support[1]):
                eq1 = Eq(R_A_x ,0)
                eq2 = Eq(R_B_x+net_x,0)
                eq3 = Eq(R_A_y + R_B_y + net_y,0)
                eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
                         R_B_x*(self._right_support[1]-self._left_support[1])+moment_A,0)
                eq5 = Eq(moment_hinge_left + R_A_y*(self._left_support[0]-self._crown_x) -\
                         T*(self._member[2]-self._crown_y),0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

            elif self._member[2]>=self._left_support[1]:
                eq1 = Eq(R_A_x ,0)
                eq2 = Eq(R_B_x+ T +net_x,0)
                eq3 = Eq(R_A_y + R_B_y + net_y,0)
                eq4 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
                         R_B_x*(self._right_support[1]-self._left_support[1])-\
                         T*(self._member[2]-self._left_support[0])+moment_A,0)
                eq5 = Eq(moment_hinge_left + R_A_y*(self._left_support[0]-self._crown_x)-\
                         T*(self._member[2]-self._crown_y),0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

            elif self._member[2]>=self._right_support[0]:
                eq1 = Eq(R_A_x,0)
                eq2 = Eq(R_B_x- T +net_x,0)
                eq3 = Eq(R_A_y + R_B_y + net_y,0)
                eq4 = Eq(moment_hinge_left+R_A_y*(self._left_support[0]-self._crown_x),0)
                eq5 = Eq(moment_A+R_B_y*(self._right_support[0]-self._left_support[0])-\
                         R_B_x*(self._right_support[1]-self._left_support[1])+\
                         T*(self._member[2]-self._left_support[1]),0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

        elif self._supports['right'] == 'roller':
            if self._member[2]>=max(self._left_support[1], self._right_support[1]):
                eq1 = Eq(R_B_x,0)
                eq2 = Eq(R_A_x+net_x,0)
                eq3 = Eq(R_A_y+R_B_y+net_y,0)
                eq4 = Eq(moment_hinge_right+R_B_y*(self._right_support[0]-self._crown_x)+\
                         T*(self._member[2]-self._crown_y),0)
                eq5 = Eq(moment_A+R_B_y*(self._right_support[0]-self._left_support[0]),0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

            elif self._member[2]>=self._left_support[1]:
                eq1 = Eq(R_B_x,0)
                eq2 = Eq(R_A_x+T+net_x,0)
                eq3 = Eq(R_A_y+R_B_y+net_y,0)
                eq4 = Eq(moment_hinge_right+R_B_y*(self._right_support[0]-self._crown_x),0)
                eq5 = Eq(moment_A-T*(self._member[2]-self._left_support[1])+\
                         R_B_y*(self._right_support[0]-self._left_support[0]),0)
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))

            elif self._member[2]>=self._right_support[1]:
                eq1 = Eq(R_B_x,0)
                eq2 = Eq(R_A_x-T+net_x,0)
                eq3 = Eq(R_A_y+R_B_y+net_y,0)
                eq4 = Eq(moment_hinge_right+R_B_y*(self._right_support[0]-self._crown_x)+\
                         T*(self._member[2]-self._crown_y),0)
                eq5 = Eq(moment_A+T*(self._member[2]-self._left_support[1])+\
                         R_B_y*(self._right_support[0]-self._left_support[0]))
                solution = solve((eq1,eq2,eq3,eq4,eq5),(R_A_x,R_A_y,R_B_x,R_B_y,T))
        else:
            eq1 = Eq(R_A_x + R_B_x + net_x,0)
            eq2 = Eq(R_A_y + R_B_y + net_y,0)
            eq3 = Eq(R_B_y*(self._right_support[0]-self._left_support[0])-\
                     R_B_x*(self._right_support[1]-self._left_support[1])+moment_A,0)
            eq4 = Eq(moment_hinge_right + R_B_y*(self._right_support[0]-self._crown_x) -\
                     R_B_x*(self._right_support[1]-self._crown_y),0)
            solution = solve((eq1,eq2,eq3,eq4),(R_A_x,R_A_y,R_B_x,R_B_y))

        for symb in self._reaction_force:
            self._reaction_force[symb] = solution[symb]

        self._bending_moment = - (self._moment_x_func.subs(x,x0) + self._moment_y_func.subs(x,x0) -\
                                  solution[R_A_y]*(x0-self._left_support[0]) +\
                                  solution[R_A_x]*(self._shape_eqn.subs({x:x0})-self._left_support[1]))

        angle  = atan(diff(self._shape_eqn,x))

        fx = (self._load_x_func+solution[R_A_x])
        fy = (self._load_y_func+solution[R_A_y])

        axial_force = fx*cos(angle) + fy*sin(angle)
        shear_force = -fx*sin(angle) + fy*cos(angle)

        self._axial_force = axial_force
        self._shear_force = shear_force

    @doctest_depends_on(modules=('numpy',))
    def draw(self):
        """
        This method returns a plot object containing the diagram of the specified arch along with the supports
        and forces applied to the structure.

        Examples
        ========

        >>> from sympy import Symbol
        >>> t = Symbol('t')
        >>> from sympy.physics.continuum_mechanics.arch import Arch
        >>> a = Arch((0,0),(40,0),crown_x=20,crown_y=12)
        >>> a.apply_load(-1,'C',8,150,angle=270)
        >>> a.apply_load(0,'D',start=20,end=40,mag=-4)
        >>> a.apply_load(-1,'E',10,t,angle=300)
        >>> p = a.draw()
        >>> p # doctest: +ELLIPSIS
        Plot object containing:
        [0]: cartesian line: 11.325 - 3*(x - 20)**2/100 for x over (0.0, 40.0)
        [1]: cartesian line: 12 - 3*(x - 20)**2/100 for x over (0.0, 40.0)
        ...
        >>> p.show()

        """
        x = Symbol('x')
        markers = []
        annotations = self._draw_loads()
        rectangles = []
        supports = self._draw_supports()
        markers+=supports

        xmax = self._right_support[0]
        xmin = self._left_support[0]
        ymin = min(self._left_support[1],self._right_support[1])
        ymax = self._crown_y

        lim = max(xmax*1.1-xmin*0.8+1, ymax*1.1-ymin*0.8+1)

        rectangles = self._draw_rectangles()

        filler = self._draw_filler()
        rectangles+=filler

        if self._member is not None:
            if(self._member[2]>=self._right_support[1]):
                markers.append(
                    {
                        'args':[[self._member[1]+0.005*lim],[self._member[2]]],
                        'marker':'o',
                        'markersize': 4,
                        'color': 'white',
                        'markerfacecolor':'none'
                    }
                )

            if(self._member[2]>=self._left_support[1]):
                markers.append(
                    {
                        'args':[[self._member[0]-0.005*lim],[self._member[2]]],
                        'marker':'o',
                        'markersize': 4,
                        'color': 'white',
                        'markerfacecolor':'none'
                    }
                )



        markers.append({
            'args':[[self._crown_x],[self._crown_y-0.005*lim]],
            'marker':'o',
            'markersize': 5,
            'color':'white',
            'markerfacecolor':'none',
        })

        if lim==xmax*1.1-xmin*0.8+1:

            sing_plot = plot(self._shape_eqn-0.015*lim,
                             self._shape_eqn,
                             (x, self._left_support[0], self._right_support[0]),
                             markers=markers,
                             show=False,
                             annotations=annotations,
                             rectangles = rectangles,
                             xlim=(xmin-0.05*lim, xmax*1.1),
                             ylim=(xmin-0.05*lim, xmax*1.1),
                             axis=False,
                             line_color='brown')

        else:
            sing_plot = plot(self._shape_eqn-0.015*lim,
                             self._shape_eqn,
                             (x, self._left_support[0], self._right_support[0]),
                             markers=markers,
                             show=False,
                             annotations=annotations,
                             rectangles = rectangles,
                             xlim=(ymin-0.05*lim, ymax*1.1),
                             ylim=(ymin-0.05*lim, ymax*1.1),
                             axis=False,
                             line_color='brown')

        return sing_plot


    def _draw_supports(self):
        support_markers = []

        xmax = self._right_support[0]
        xmin = self._left_support[0]
        ymin = min(self._left_support[1],self._right_support[1])
        ymax = self._crown_y

        if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
            max_diff = 1.1*xmax-0.8*xmin
        else:
            max_diff = 1.1*ymax-0.8*ymin

        if self._supports['left']=='roller':
            support_markers.append(
                {
                    'args':[
                        [self._left_support[0]],
                        [self._left_support[1]-0.02*max_diff]
                    ],
                    'marker':'o',
                    'markersize':11,
                    'color':'black',
                    'markerfacecolor':'none'
                }
            )
        else:
            support_markers.append(
                {
                    'args':[
                        [self._left_support[0]],
                        [self._left_support[1]-0.007*max_diff]
                    ],
                    'marker':6,
                    'markersize':15,
                    'color':'black',
                    'markerfacecolor':'none'
                }
            )

        if self._supports['right']=='roller':
            support_markers.append(
                {
                    'args':[
                        [self._right_support[0]],
                        [self._right_support[1]-0.02*max_diff]
                    ],
                    'marker':'o',
                    'markersize':11,
                    'color':'black',
                    'markerfacecolor':'none'
                }
            )
        else:
            support_markers.append(
                {
                    'args':[
                        [self._right_support[0]],
                        [self._right_support[1]-0.007*max_diff]
                    ],
                    'marker':6,
                    'markersize':15,
                    'color':'black',
                    'markerfacecolor':'none'
                }
            )

        support_markers.append(
            {
                'args':[
                    [self._right_support[0]],
                    [self._right_support[1]-0.036*max_diff]
                ],
                'marker':'_',
                'markersize':15,
                'color':'black',
                'markerfacecolor':'none'
            }
        )

        support_markers.append(
            {
                'args':[
                    [self._left_support[0]],
                    [self._left_support[1]-0.036*max_diff]
                ],
                'marker':'_',
                'markersize':15,
                'color':'black',
                'markerfacecolor':'none'
            }
        )

        return support_markers

    def _draw_rectangles(self):
        member = []

        xmax = self._right_support[0]
        xmin = self._left_support[0]
        ymin = min(self._left_support[1],self._right_support[1])
        ymax = self._crown_y

        if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
            max_diff = 1.1*xmax-0.8*xmin
        else:
            max_diff = 1.1*ymax-0.8*ymin

        if self._member is not None:
            if self._member[2]>= max(self._left_support[1],self._right_support[1]):
                member.append(
                    {
                        'xy':(self._member[0],self._member[2]-0.005*max_diff),
                        'width':self._member[1]-self._member[0],
                        'height': 0.01*max_diff,
                        'angle': 0,
                        'color':'brown',
                    }
                )

            elif self._member[2]>=self._left_support[1]:
                member.append(
                    {
                        'xy':(self._member[0],self._member[2]-0.005*max_diff),
                        'width':self._right_support[0]-self._member[0],
                        'height': 0.01*max_diff,
                        'angle': 0,
                        'color':'brown',
                    }
                )

            else:
                member.append(
                    {
                        'xy':(self._member[1],self._member[2]-0.005*max_diff),
                        'width':abs(self._left_support[0]-self._member[1]),
                        'height': 0.01*max_diff,
                        'angle': 180,
                        'color':'brown',
                    }
                )

        if self._distributed_loads:
            for loads in self._distributed_loads:

                start = self._distributed_loads[loads]['start']
                end = self._distributed_loads[loads]['end']

                member.append(
                    {
                        'xy':(start,self._crown_y+max_diff*0.15),
                        'width': (end-start),
                        'height': max_diff*0.01,
                        'color': 'orange'
                    }
                )


        return member

    def _draw_loads(self):
        load_annotations = []

        xmax = self._right_support[0]
        xmin = self._left_support[0]
        ymin = min(self._left_support[1],self._right_support[1])
        ymax = self._crown_y

        if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
            max_diff = 1.1*xmax-0.8*xmin
        else:
            max_diff = 1.1*ymax-0.8*ymin

        for load in self._conc_loads:
            x = self._conc_loads[load]['x']
            y = self._conc_loads[load]['y']
            angle = self._conc_loads[load]['angle']
            mag = self._conc_loads[load]['mag']
            load_annotations.append(
                {
                    'text':'',
                    'xy':(
                        x+cos(rad(angle))*max_diff*0.08,
                        y+sin(rad(angle))*max_diff*0.08
                    ),
                    'xytext':(x,y),
                    'fontsize':10,
                    'fontweight': 'bold',
                    'arrowprops':{'width':1.5, 'headlength':5, 'headwidth':5, 'facecolor':'blue','edgecolor':'blue'}
                }
            )
            load_annotations.append(
                {
                    'text':f'{load}: {mag} N',
                    'fontsize':10,
                    'fontweight': 'bold',
                    'xy': (x+cos(rad(angle))*max_diff*0.12,y+sin(rad(angle))*max_diff*0.12)
                }
            )

        for load in self._distributed_loads:
            start = self._distributed_loads[load]['start']
            end = self._distributed_loads[load]['end']
            mag = self._distributed_loads[load]['f_y']
            x_points = numpy.arange(start,end,(end-start)/(max_diff*0.25))
            x_points = numpy.append(x_points,end)
            for point in x_points:
                if(mag<0):
                    load_annotations.append(
                        {
                            'text':'',
                            'xy':(point,self._crown_y+max_diff*0.05),
                            'xytext': (point,self._crown_y+max_diff*0.15),
                            'arrowprops':{'width':1.5, 'headlength':5, 'headwidth':5, 'facecolor':'orange','edgecolor':'orange'}
                        }
                    )
                else:
                    load_annotations.append(
                        {
                            'text':'',
                            'xy':(point,self._crown_y+max_diff*0.2),
                            'xytext': (point,self._crown_y+max_diff*0.15),
                            'arrowprops':{'width':1.5, 'headlength':5, 'headwidth':5, 'facecolor':'orange','edgecolor':'orange'}
                        }
                    )
            if(mag<0):
                load_annotations.append(
                    {
                        'text':f'{load}: {abs(mag)} N/m',
                        'fontsize':10,
                        'fontweight': 'bold',
                        'xy':((start+end)/2,self._crown_y+max_diff*0.175)
                    }
                )
            else:
                load_annotations.append(
                    {
                        'text':f'{load}: {abs(mag)} N/m',
                        'fontsize':10,
                        'fontweight': 'bold',
                        'xy':((start+end)/2,self._crown_y+max_diff*0.125)
                    }
                )
        return load_annotations

    def _draw_filler(self):
        x = Symbol('x')
        filler = []
        xmax = self._right_support[0]
        xmin = self._left_support[0]
        ymin = min(self._left_support[1],self._right_support[1])
        ymax = self._crown_y

        if abs(1.1*xmax-0.8*xmin)>abs(1.1*ymax-0.8*ymin):
            max_diff = 1.1*xmax-0.8*xmin
        else:
            max_diff = 1.1*ymax-0.8*ymin

        x_points = numpy.arange(self._left_support[0],self._right_support[0],(self._right_support[0]-self._left_support[0])/(max_diff*max_diff))

        for point in x_points:
            filler.append(
                    {
                        'xy':(point,self._shape_eqn.subs(x,point)-max_diff*0.015),
                        'width': (self._right_support[0]-self._left_support[0])/(max_diff*max_diff),
                        'height': max_diff*0.015,
                        'color': 'brown'
                    }
            )

        return filler
