@@ -2999,7 +2999,7 @@ def SHarmcal(SytSym,SHFln,psi,gam):
29992999 SHVal += (SHFln [term ][0 ]* Ksl )
30003000 return SHVal
30013001
3002- def KslCalc (trm ,psi ,gam ):
3002+ def KslCalc (trm ,psi ,gam , ifDerv = False ):
30033003 '''Compute one angular part term in spherical harmonics
30043004
30053005 :param str trm:sp. harm term name in the form of 'C(l,m)' or 'C(l,m)c' for cubic
@@ -3008,11 +3008,22 @@ def KslCalc(trm,psi,gam):
30083008
30093009 :returns array Ksl: spherical harmonics angular part for psi,gam pairs
30103010 '''
3011+ delt = 0.0001
30113012 l ,m = eval (trm .strip ('C' ).strip ('c' ))
30123013 if 'c' in trm :
3013- return CubicSHarm (l ,m ,psi ,gam )
3014+ if ifDerv :
3015+ dYdp = (CubicSHarm (l ,m ,psi + delt ,gam )- CubicSHarm (l ,m ,psi - delt ,gam ))/ (2.0 * delt )
3016+ dYdg = (CubicSHarm (l ,m ,psi ,gam + delt )- CubicSHarm (l ,m ,psi ,gam - delt ))/ (2.0 * delt )
3017+ return CubicSHarm (l ,m ,psi ,gam ),dYdp ,dYdg
3018+ else :
3019+ return CubicSHarm (l ,m ,psi ,gam )
30143020 else :
3015- return SphHarmAng (l ,m ,1.0 ,psi ,gam )
3021+ if ifDerv :
3022+ dYdp = (SphHarmAng (l ,m ,1.0 ,psi + delt ,gam )- SphHarmAng (l ,m ,1.0 ,psi - delt ,gam ))/ (2.0 * delt )
3023+ dYdg = (SphHarmAng (l ,m ,1.0 ,psi ,gam + delt )- SphHarmAng (l ,m ,1.0 ,psi ,gam - delt ))/ (2.0 * delt )
3024+ return SphHarmAng (l ,m ,1.0 ,psi ,gam ),dYdp ,dYdg
3025+ else :
3026+ return SphHarmAng (l ,m ,1.0 ,psi ,gam )
30163027
30173028def SphHarmAng (L ,M ,P ,Th ,Ph ):
30183029 ''' Compute spherical harmonics values using scipy.special.sph_harm
0 commit comments