#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Package DLMcalcs contains functions to load the influence coeffient matrix
calculation c function and other necessary calculations for the Doublet Lattice
Method.

This code is part of the SDPMflut Python distribution.
Copyright (C) 2024 Grigorios Dimitriadis

This program is free software: you can redistribute it and/or modify it 
under the terms of the GNU General Public License as published by the 
Free Software Foundation, either version 3 of the License, or (at your 
option) any later version.
 
This program is distributed in the hope that it will be useful, but 
WITHOUT ANY WARRANTY; without even the implied warranty of 
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General 
Public License for more details.

You should have received a copy of the GNU General Public License along 
with this program. If not, see <https://www.gnu.org/licenses/>.
"""

import ctypes
import numpy as np
import platform
from numpy.ctypeslib import ndpointer
from scipy import linalg
import numpy.matlib

def Cfuncs(install_dir):
    # Loads the c function used by the nonplanar DLM included in SDPMflut. 
    # install_dir: Path to ./Common directory
    # dlminf_nonplanarso: function to calculate the influence coefficient matrix
    
    # Load c function to calculate steady SDPM influence coefficients 
    if platform.system() == 'Windows':
        so_file = install_dir+r"\dlminf_nonplanarso.so" 
    else:
        so_file = install_dir+r"dlminf_nonplanarso.so" 
    # End if
    lib = ctypes.cdll.LoadLibrary(so_file)
    dlminf_nonplanar = lib.dlminf_nonplanar
    dlminf_nonplanar.restype = None
    # Define arguments of sdpminf
    dlminf_nonplanar.argtypes = [ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),               
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),               
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),               
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),               
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),               
                    ndpointer(ctypes.c_double, flags="C_CONTIGUOUS"),               
    		        ctypes.c_int,ctypes.c_int,ctypes.c_int,ctypes.c_int]
    
    return dlminf_nonplanar

def DLMdtypes():
    # Defines all data types used by the DLM method in SDPMflut. These are 
    # numpy.dtype variables for numpy arrays.
    # tp_trap: Data type for defining trapezoidal sections of wings
    # tp_body: Data type for defining SDPM grids for bodies
    # tp_allbodies: Data type for concatenating grid and result information for 
    #               all bodies in a flow
    
    # Define trapezoidal section data type
    tp_trap=np.dtype({'names':('rootchord', 'xledist','span', 'taper','sweepLE',
                               'roottwist','tiptwist','twistcent','dihedral','rootairfoil',
                               'rootairfoilparams','tipairfoil','tipairfoilparams'),
              'formats':('f8', 'f8', 'f8', 'f8','f8', 'f8', 'f8', 'f8','f8','U25','(2,)f8','U25','(2,)f8')})
    # Define body data type
    tp_body=np.dtype({'names':('Xp0', 'Yp0', 'Zp0', 
                          'Xc0','Yc0','Zc0',
                          'Xc0all','Yc0all','Zc0all', 
                          'Xm','Ym','Zm','s',
                          'Xd','Yd','Zd','s0', 
                          'Deltax','Deltay',
                          'm','n','c0','yc', 
                          'AR','S','name','Schar', 
                          'Phi_x','Phi_y','Phi_z', 
                          'Phi_phi','Phi_theta','Phi_psi', 
                          'Phi_xall','Phi_yall','Phi_zall', 
                          'Phi_phiall','Phi_thetaall','Phi_psiall', 
                          'Xp', 'Yp', 'Zp', 
                          'Xc','Yc', 'Zc',
                          'Xcall','Ycall','Zcall',
                          'nxall','nyall','nzall',
                          'nx','ny','nz',
                          'Zccamb_diff','Zccamb_diffall','gammas',
                          's0all','sall','Deltacp0', 
                          'Fx0','Fy0','Fz0','bchar',
                          'Mx0','My0','Mz0','b'),    
              'formats':(object,object,object,                          #'Xp0', 'Yp0', 'Zp0'
                        object,object,object,                           #'Xc0','Yc0','Zc0'
                        object,object,object,                           # 'Xc0all','Yc0all','Zc0all'
                        object,object,object,object,                    # 'Xm','Ym','Zm','s'
                        object,object,object,object,                    # 'Xd','Yd','Zd','s0'
                        object,object,                                  # 'Deltax','Deltay'
                        np.int64,np.int64,np.float64,object,            # 'm','n','c0','yc'
                        np.float64,np.float64,'S25',np.float64,         # 'AR','S','name','Schar'
                        object,object,object,                           # 'Phi_x','Phi_y','Phi_z'
                        object,object,object,                           # 'Phi_phi','Phi_theta','Phi_psi'
                        object,object,object,                           # 'Phi_xall','Phi_yall','Phi_zall'
                        object,object,object,                           # 'Phi_phiall','Phi_thetaall','Phi_psiall'
                        object,object,object,                           # 'Xp','Yp','Zp'
                        object,object,object,                           # 'Xc','Yc','Zc',
                        object,object,object,                           # 'Xcall','Ycall','Zcall'
                        object,object,object,                           # 'nxall','nyall','nzall'                        
                        object,object,object,                           # 'nx','ny','nz'                        
                        object,object,object,                           # 'Zccamb_diff','Zccamb_diffall','gammas'
                        object,object,object,                           # s0all','sall','Deltacp0'
                        object,object,object,np.float64,                # 'Fx0','Fy0','Fz0','bchar'
                        object,object,object,np.float64)})              # 'Mx0','My0','Mz0','b'    
    # Define allbodies data type
    tp_allbodies=np.dtype({'names':('nbody', 'bodypanels', 'bodypanelsn', 
                          'allpanelsn','allpanels',
                          'inds','bchar', 
                          'Xc0all','Yc0all','Zc0all',
                          'nxall','nyall','nzall',
                          's0all','Zccamb_diffall', 
                          'Xcall','Ycall','Zcall', 
                          'Phi_xall','Phi_yall','Phi_zall', 
                          'Phi_phiall','Phi_thetaall','Phi_psiall', 
                          'Deltacp0','Fx0','Fy0','Fz0',
                          'Mx0','My0','Mz0'),    
              'formats':(np.int64,object,object,                    # 'nbody', 'bodypanels', 'bodypanelsn'
                        np.int64,np.int64,                          # 'allpanelsn','allpanels'
                        object,np.float64,                          # 'inds','bchar' 
                        object,object,object,                       # 'Xc0all','Yc0all','Zc0all',
                        object,object,object,                       # 'nxall','nyall','nzall',
                        object,object,                              # 's0all','Zccamb_diffall'
                        object,object,object,                       # 'Xcall','Ycall','Zcall'
                        object,object,object,                       # 'Phi_xall','Phi_yall','Phi_zall'
                        object,object,object,                       # 'Phi_phiall','Phi_thetaall','Phi_psiall'
                        object,object,object,object,                # 'Deltacp0','Fx0','Fy0','Fz0' 
                        object,object,object)})                     # 'Mx0','My0','Mz0'         
    return tp_trap, tp_body, tp_allbodies

def allbodyindexDLM(body,bchar):
    # Calculates the numbers and indices of panels, spanwise panels and wake
    # panels in all bodies stored in struct array body for the DLM.
    # body: struct array containing the geometry of all the bodies in the flow.
    #       Each element of body describes a different wing, fairing or 
    #       fuselage
    # allbodies: struct array containing concatenated information for all the
    #       bodies
   
    # Total number of bodies
    nbody=len(body)
    # Panel numbers and indexing
    bodypanels=body['m']*body['n']   # Number of panels in each body
    bodypanelsn=body['n']               # Number of spanwise panels in each body
    allpanelsn=np.sum(body['n'])        # Total number of spanwise body panels
    allpanels=np.sum(bodypanels)        # Total number of body panels
    # Indices for assigning body-on-body influence coefficient matrices to assembled arrays
    inds=np.zeros(nbody+1,dtype=int)
    inds[1:nbody+1]=np.cumsum(bodypanels) 
    
    Zccamb_diffall=np.zeros((allpanels,1)) 
    nxall=np.zeros((allpanels,1)) 
    nyall=np.zeros((allpanels,1)) 
    nzall=np.zeros((allpanels,1)) 
    Xcall=np.zeros((allpanels,1)) 
    Ycall=np.zeros((allpanels,1)) 
    Zcall=np.zeros((allpanels,1)) 
    Xc0all=np.zeros((allpanels,1)) 
    Yc0all=np.zeros((allpanels,1)) 
    Zc0all=np.zeros((allpanels,1)) 
    s0all=np.zeros((allpanels,1)) 
    # Cycle through all the bodies
    for i in range(0,nbody):
        # Camber distribution
        Zccamb_diffall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Zccamb_diffall'][i])
        # Panel normal vector components in normalized coordinates
        nxall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['nxall'][i])
        nyall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['nyall'][i])
        nzall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['nzall'][i])
        # Panel control point coordinates in normalized coordinates
        Xcall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Xcall'][i])
        Ycall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Ycall'][i])
        Zcall[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Zcall'][i])
        # Panel control point coordinates in physical coordinates
        Xc0all[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Xc0all'][i])
        Yc0all[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Yc0all'][i])
        Zc0all[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['Zc0all'][i])
        # Panel areas in physical coordinates
        s0all[inds[i]:inds[i+1],0]=np.ndarray.flatten(body['s0all'][i])
    # End for    

    # Obtain tp_allbodies data type
    _, _, tp_allbodies=DLMdtypes()
    # Initialize struct array allbodies 
    allbodies=np.zeros(1,dtype=tp_allbodies)
    # Assign panel numbers and indices. The rest of the information in 
    # allbodies will be assigned by other functions.
    allbodies[0]['nbody']=nbody
    allbodies[0]['bodypanels']=bodypanels
    allbodies[0]['bodypanelsn']=bodypanelsn
    allbodies[0]['allpanelsn']=allpanelsn
    allbodies[0]['allpanels']=allpanels
    allbodies[0]['inds']=inds
    allbodies[0]['nxall']=nxall
    allbodies[0]['nyall']=nyall
    allbodies[0]['nzall']=nzall
    allbodies[0]['Xcall']=Xcall
    allbodies[0]['Ycall']=Ycall
    allbodies[0]['Zcall']=Zcall
    allbodies[0]['Xc0all']=Xc0all
    allbodies[0]['Yc0all']=Yc0all    
    allbodies[0]['Zc0all']=Zc0all    
    allbodies[0]['s0all']=s0all    
    allbodies[0]['Zccamb_diffall']=Zccamb_diffall    
    allbodies[0]['bchar']=bchar    
    
    return allbodies    

def infcoef(body,allbodies,install_dir,k,Mach):
    # Calculate doublet lattice influence coefficient matrix.
    # body: struct array containing the geometry of all the bodies in the flow.
    #       Each element of body describes a different wing, fairing or 
    #       fuselage
    # allbodies: struct array containing concatenated information for all the
    #       bodies   
    # install_dir: Path to ./Common directory
    # k: reduced frequency value
    # Mach: Mach number
    # A: Complex influence coefficient matrix
    
    # Load c function to calculate the steady SDPM influence coefficient
    # matrices
    dlminf_nonplanar=Cfuncs(install_dir)
    # Set reduced frequency and Mach values
    params=np.array([k, Mach])                
    
    # Obtain number of bodies
    nbody=allbodies['nbody'][0]
    # Obtain total number of panels in all bodies
    allpanels=allbodies['allpanels'][0]
    # Obtain start and end indices of panels of each body
    inds=allbodies['inds'][0]
    # Initialize influence coefficient matrix
    A=np.zeros((allpanels,allpanels),dtype=complex)
    
    # Cycle through the bodies
    for i in range(0,nbody):
        for j in range(0,nbody):
            # Initialize body source influence coefficient matrix
            Areal=np.zeros((body['m'][i]*body['n'][i],body['m'][j]*body['n'][j]))
            # Initialize body doublet influence coefficient matrix
            Aimag=np.zeros((body['m'][i]*body['n'][i],body['m'][j]*body['n'][j]))
            # Calculate real and imaginary parts of influence coefficient matrix
            dlminf_nonplanar(body['Xc'][i],body['Yc'][i],body['Zc'][i],body['Xd'][j],body['Yd'][j],body['Zd'][j],
                    body['Xm'][j],body['Ym'][j],body['Zm'][j],body['gammas'][i],body['gammas'][j],
                    body['Deltax'][j],body['Deltay'][j],params,Areal,Aimag,body['m'][i],body['n'][i],body['m'][j],body['n'][j])
            # Assemble influence coefficient matrix
            A[inds[i]:inds[i+1],inds[j]:inds[j+1]]=Areal+1j*Aimag
        # End for
    # End for
    return A

def aeroforcesDLM(Deltacp,nxall,nyall,nzall,s0all,Xc0all,Yc0all,Zc0all,xf0,yf0,zf0):
    # Calculates aerodynamic loads acting on all panels. The moments
    # are taken around point xf0,yf0,zf0
    # Deltacp: Pressure jump across all panels
    # nxall, nyall, nzall: x, y and z components of unit vectors normal to all
    #      the panels.
    # s0all: areas of all the panels in physical coordinates.
    # Xc0all, Yc0all, Zc0all: x, y and z coordinates of the control points of
    #       all the panels in physical coordinates.
    # xf0, yf0, zf0: x, y and z coordinates of the point around which moments
    #       are to be taken in physical coordinates.
    # Fx, Fy, Fz: x, y and z components of aerodynamic force
    # Mx, My, Mz: aerodynamic moments around the x, y and z axes.
    
    Fx=Deltacp*nxall*s0all   # Force in x direction
    Fy=Deltacp*nyall*s0all   # Force in x direction
    Fz=Deltacp*nzall*s0all   # Force in x direction
    Mx=(Yc0all-yf0)*Fz-(Zc0all-zf0)*Fy   # Rolling moment
    My=-(Xc0all-xf0)*Fz+(Zc0all-zf0)*Fx  # Pitching moment
    Mz=(Xc0all-xf0)*Fy-(Yc0all-yf0)*Fx   # Yawing moment  
    
    return Fx,Fy,Fz,Mx,My,Mz

def steadysolveDLM(body,allbodies,Mach,alpha0,beta0,xf0,yf0,zf0,install_dir):
    # Calculates the steady aerodynamic pressure and loads acting on the
    # panel control points of the bodies described in struct array body, using
    # the Doublet Lattice Method.
    # body: struct array containing the geometry of all the bodies in the flow.
    #       Each element of body describes a different wing, fairing or 
    #       fuselage
    # allbodies: struct array containing concatenated information for all the
    #       bodies
    # Mach: Free stream Mach number
    # alpha0: Free stream angle of attack
    # beta0: Free stream sideslip angle
    # xf0,yf0,zf0: Chordwise and spanwise cartesian coordinates of centre of 
    # gravity or rotation
    # install_dir: Path to ./Common directory
    # Cx0, Cy0, Cz0: Aerodynamic force coefficients in x, y and z directions.
    # Cl0, Cm0, Cz0: Aerodynamic moment coefficients around x, y and z axes.
    
    # Set reduced frequency to zero for steady flow
    k=0.0
    # Calculate influence coefficient matrix
    Q=infcoef(body,allbodies,install_dir,k,Mach)
    
    # Determine number of angles of attack to calculate
    nalpha=alpha0.size
    # Initialize result arrays
    Cx0=np.zeros(nalpha)
    Cy0=np.zeros(nalpha)
    Cz0=np.zeros(nalpha)
    Cl0=np.zeros(nalpha)
    Cm0=np.zeros(nalpha)
    Cn0=np.zeros(nalpha)
    # Cycle through the angles of attack
    for ialpha in range(0,nalpha):
        
        # Calculate non-dimensional free stream velocity components
        barUinf=np.cos(alpha0[ialpha])*np.cos(beta0)
        barVinf=-np.sin(beta0)
        barWinf=np.sin(alpha0[ialpha])*np.cos(beta0)
        # Calculate upwash
        w0=(barVinf-allbodies[0]['Zccamb_diffall'])*allbodies[0]['nyall'] \
            +(barWinf-allbodies[0]['Zccamb_diffall'])*allbodies[0]['nzall']
        # Calculate steady pressure jump coefficient on control points
        Deltacp0=np.real(linalg.solve(-Q,w0)) # The result will be real since k=0
        # Calculate aerodynamic loads on control points
        Fx0,Fy0,Fz0,Mx0,My0,Mz0=aeroforcesDLM(Deltacp0,allbodies[0]['nxall'],allbodies[0]['nyall'],allbodies[0]['nzall'], \
               allbodies[0]['s0all'],allbodies[0]['Xc0all'],allbodies[0]['Yc0all'],allbodies[0]['Zc0all'],xf0,yf0,zf0)

        # Assign concatenated steady pressure coefficient and aerodynamic loads to
        # struct array allbodies
        allbodies[0]['Deltacp0']=Deltacp0
        allbodies[0]['Fx0']=Fz0
        allbodies[0]['Fy0']=Fz0
        allbodies[0]['Fz0']=Fz0
        allbodies[0]['Mx0']=Mx0
        allbodies[0]['My0']=My0
        allbodies[0]['Mz0']=My0
        # Assign steady pressure coefficient and aerodynamic loads to respective 
        # elements of struct array body
        for i in range(0,allbodies['nbody'][0]):
            body['Deltacp0'][i]=np.reshape(Deltacp0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
            body['Fx0'][i]=np.reshape(Fx0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
            body['Fy0'][i]=np.reshape(Fy0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
            body['Fz0'][i]=np.reshape(Fz0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
            body['Mx0'][i]=np.reshape(Mx0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
            body['My0'][i]=np.reshape(My0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
            body['Mz0'][i]=np.reshape(Mz0[allbodies['inds'][0][i]:allbodies['inds'][0][i+1],0],(body['m'][i],body['n'][i]),order='C')
        # End if 
    
        # Sum up aerodynamic loads
        Cx0[ialpha]=np.sum(Fx0)/(2.0*body[0]['S'])
        Cy0[ialpha]=np.sum(Fy0)/(2.0*body[0]['S'])
        Cz0[ialpha]=np.sum(Fz0)/(2.0*body[0]['S'])
        Cl0[ialpha]=np.sum(Mx0)/(2.0*body[0]['S'])/body[0]['b']
        Cm0[ialpha]=np.sum(My0)/(2.0*body[0]['S'])/allbodies[0]['bchar']
        Cn0[ialpha]=np.sum(Mz0)/(2.0*body[0]['S'])/body[0]['b']
    # End for
    
    return body,allbodies,barUinf,barVinf,barWinf,Cx0,Cy0,Cz0,Cl0,Cm0,Cn0

def unsteadysolve_flexDLM(body,allbodies,k,Mach,install_dir):
    # Calculates the unsteady aerodynamic pressure and loads acting on the
    # panel control points of the bodies described in struct array body, using 
    # the Doublet Lattice Method. The motion of the bodies is flexible and 
    # parallel to pre-defined mode shapes.
    # body: struct array containing the geometry of all the bodies in the flow.
    #       Each element of body describes a different wing, fairing or 
    #       fuselage
    # allbodies: struct array containing concatenated information for all the
    #       bodies
    # k: Reduced frequency
    # Mach: Free stream Mach number
    # install_dir: Path to ./Common directory
    # Deltacp1: Total pressure coefficient per generalized coordinate acting on the 
    #       control points of all bodies in the flow
    # Deltacp_0: Stiffness pressure coefficient per generalized coordinate acting on  
    #       the control points of all bodies in the flow
    # Deltacp_1: Damping pressure coefficient per generalized coordinate acting on  
    #       the control points of all bodies in the flow
    
    # Calculate influence coefficient matrix
    Q=infcoef(body,allbodies,install_dir,k,Mach)
    
    # Calculate pressure derivatives with respect to modes of vibration
    # Displacements    
    Deltacptheta=-linalg.solve(Q,allbodies[0]['Phi_thetaall']*allbodies[0]['nzall'])
    Deltacppsi=-linalg.solve(Q,-allbodies[0]['Phi_psiall']*allbodies[0]['nyall'])
    Deltacpydot=-linalg.solve(Q,-1/allbodies[0]['bchar']*allbodies[0]['Phi_yall']*allbodies[0]['nyall'])
    Deltacpzdot=-linalg.solve(Q,-1/allbodies[0]['bchar']*allbodies[0]['Phi_zall']*allbodies[0]['nzall'])

    # Calculate oscillatory pressure components
    Deltacp_0=Deltacptheta+Deltacppsi    # Displacement component
    Deltacp_1=Deltacpydot+Deltacpzdot  # Velocity component
    # Calculate total oscillatory pressure
    Deltacp1=Deltacp_0+(1j*k)*Deltacp_1

    return Deltacp1,Deltacp_0,Deltacp_1
    
def unsteadysolve_pitchplungeDLM(body,allbodies,k,Mach,xf0,yf0,install_dir):
    # Calculates the unsteady aerodynamic pressure and loads acting on the
    # panel control points of the bodies described in struct array body, using 
    # the Doublet Lattice Method. The motion of the bodies is rigid pitch and 
    # plunge.
    # body: struct array containing the geometry of all the bodies in the flow.
    #       Each element of body describes a different wing, fairing or 
    #       fuselage
    # allbodies: struct array containing concatenated information for all the
    #       bodies
    # k: Reduced frequency
    # Mach: Free stream Mach number
    # xf,yf,zf: Position of pitch axis in Prandtl-Glauert coordinates
    # install_dir: Path to ./Common directory
    # cpalpha: Pressure derivative with respect to pitch
    # cphdot: Pressure derivative with respect to plunge velocity
    # cpalphadot: Pressure derivative with respect to pitch velocity
    # cph2dot: Pressure derivative with respect to plunge acceleration
    # cpalpha2dot: Pressure derivative with respect to pitch acceleration

    # Calculate influence coefficient matrix
    Q=infcoef(body,allbodies,install_dir,k,Mach)
 
    # Calculate pressure jump derivatives
    Deltacphdot=-linalg.solve(Q,(1.0/allbodies[0]['bchar']*np.ones((allbodies[0]['allpanels'],1))))
    Deltacpalpha=-linalg.solve(Q,np.ones((allbodies[0]['allpanels'],1)))
    Deltacpalphadot=-linalg.solve(Q,(allbodies[0]['Xcall']-xf0/allbodies[0]['bchar']))
    
    return Deltacphdot,Deltacpalpha,Deltacpalphadot
        
