dipole_integral

Dipole integral sample.

''' Dipole integral sample.
'''
from irfpy.util import dipole
import numpy as np

def main():
    '''Main script'''

    ### Dipole traced. L=10.

    lvalue = 10
#    ndiv = 9
    ndiv = 90
    ndiv2 = ndiv * 2 + 1
    r, tr = dipole.fieldline_rt(lvalue, ndiv=ndiv2)
    t = np.rad2deg(tr)

    import matplotlib.pyplot as plt
    plt.subplot(221)
    plt.plot(t, r, 'o')
    plt.xlabel('Angle')
    plt.ylabel('Distance')

    x = r * np.sin(tr)
    z = r * np.cos(tr)

    plt.subplot(222)
    plt.plot(x, z, 'o')
    plt.xlabel('x')
    plt.ylabel('z')

    ### Numerical integral

    # Calculate distance in between
    dx = x[1:] - x[:-1]
    dz = z[1:] - z[:-1]
    dl = np.sqrt(dx ** 2 + dz ** 2)

    l = np.cumsum(dl)

    plt.subplot(223)
    plt.plot(t[1:], l, 'o')

    ### My expression
    ## value at 0 deg
    dm0 = distance_mother(0, lvalue=lvalue)
    print(dm0)

    distlist = []
    for dmtr in range(ndiv2):
        dm = distance_mother(tr[dmtr], lvalue=lvalue)
        distlist.append(dm - dm0)

    plt.plot(t, distlist, 'rx')
        
    plt.savefig('dipole_integral.png')

    
def distance_mother(theta, lvalue=1.0):
    ''' theta in radian.
    '''
    t = np.cos(theta)
    val = - lvalue * ( 1 / 6. * (3 * np.sqrt(3 * (t ** 2) + 1) * t + np.sqrt(3) * np.arcsinh(np.sqrt(3) * t)))
    
    return val
    

    



if __name__ == "__main__":
    main()