<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Index
Weak-field approximation for Stokes V
Strong-field approximation for Stokes V
Longitudinal magnetograph
Center-of-gravity
Unresolved fields - inclination angle
Unresolved fields - line ratio
Weak-field approximation Stokes V
<a id='weak_field_V'></a>
Let us use the linear relation between the circular polarization profile and the derivative of the intensity profile. This amounts to solving a linear fit problem, that can be analytically obtained
Step1:  Strong-field approximation Stokes V
<a id='strong_field_V'></a>
In the strong field approximation, the separation of the components is just proportional to $\Delta \lambda_B$. So we only need to calculate this separation and get the field.
Step2:  Longitudinal magnetograph
<a id='longMagnetograph'></a>
Step3:  Center-of-gravity
<a id='cog'></a>
Step4:  Influence of filling factor on inclination of magnetic field
<a id='inclinationFilling'></a>
Step5:  Line ratio technique
<a id='lineRatio'></a>
<ASSISTANT_TASK:>
Python Code:
lambda0 = 6301.5
JUp = 1.0
JLow = 1.0
gUp = 2.5
gLow = 0.0
lambdaStart = 6300.8
lambdaStep = 0.01
nLambda = 150
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep
lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])
s = pymilne.milne(nLambda, lineInfo)
BField = 100.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.080
kl = 6.75
modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes = s.synth(modelSingle,mu)
cont = B0+B1
f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
ax[0].plot(wavelength, stokes[0,:] / cont)
ax[1].plot(wavelength, stokes[3,:] / cont)
ax[0].set_xlabel('Wavelength [$\AA$]')
ax[1].set_xlabel('Wavelength [$\AA$]')
ax[0].set_ylabel('I/I$_c$')
ax[1].set_ylabel('V/I$_c$')
ax[0].ticklabel_format(useOffset=False)
ax[1].ticklabel_format(useOffset=False)
pl.tight_layout()
geff = 0.5*(gUp+gLow) + 0.25*(gUp-gLow)*(JUp*(JUp+1)-JLow*(JLow+1))
diffI = np.gradient(stokes[0,:], lambdaStep)
alpha = -4.6686e-13*6301.5**2*geff
print "Estimated field = {0} G".format(np.sum(stokes[3,:] * diffI) / np.sum(diffI**2) / alpha)
Explanation: Index
Weak-field approximation for Stokes V
Strong-field approximation for Stokes V
Longitudinal magnetograph
Center-of-gravity
Unresolved fields - inclination angle
Unresolved fields - line ratio
Weak-field approximation Stokes V
<a id='weak_field_V'></a>
Let us use the linear relation between the circular polarization profile and the derivative of the intensity profile. This amounts to solving a linear fit problem, that can be analytically obtained:
$$
V = \alpha \frac{dI}{d\lambda}
$$
This is equivalent to doing the linear fit $y=ax$, which can be trivially solved by writing the $\chi^2$ merit function, differentiating it and equating to zero:
$$
\chi^2 = \sum_i (y_i-ax_i)^2
$$
We solve
$$
\frac{d \chi^2}{da} = -2\sum_i (y_i-a x_i)x_i=0,
$$
with the solution
$$
a = \frac{\sum_i y_i x_i}{\sum_i x_i^2}
$$
End of explanation
BField = 10500.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.080
kl = 6.75
modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes = s.synth(modelSingle,mu)
cont = B0+B1
f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
ax[0].plot(wavelength, stokes[0,:] / cont)
ax[1].plot(wavelength, stokes[3,:] / cont)
ax[0].set_xlabel('Wavelength [$\AA$]')
ax[1].set_xlabel('Wavelength [$\AA$]')
ax[0].set_ylabel('I/I$_c$')
ax[1].set_ylabel('V/I$_c$')
ax[0].ticklabel_format(useOffset=False)
ax[1].ticklabel_format(useOffset=False)
pl.tight_layout()
maxLoc = np.argmax(stokes[3,:])
deltaL = lambda0 - wavelength[maxLoc]
print "Estimated field = {0} G".format(deltaL / (4.6686e-13*6301.**2*gUp))
Explanation: Strong-field approximation Stokes V
<a id='strong_field_V'></a>
In the strong field approximation, the separation of the components is just proportional to $\Delta \lambda_B$. So we only need to calculate this separation and get the field.
End of explanation
lambda0 = 6302.5
JUp = 1.0
JLow = 1.0
gUp = 2.5
gLow = 0.0
lambdaStart = 6301.8
lambdaStep = 0.03
nLambda = 50
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep
lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])
s = pymilne.milne(nLambda, lineInfo)
BField = 100.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.040
kl = 6.75
pV = np.exp(-(wavelength-lambda0-60e-3)**2 / 75e-3**2)
pQ = np.exp(-(wavelength-lambda0-120e-3)**2 / 75e-3**2)
nB = 100
B = np.linspace(0.0,3500.0,nB)
Bpar = np.zeros((2,4,nB))
Bperp = np.zeros((4,nB))
SV = np.zeros((2,4,nB))
SQ = np.zeros((4,nB))
pal = sn.color_palette()
thetas = [0, 30, 60, 85]
for j in range(4):
	for i in range(nB):
		modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
		stokes = s.synth(modelSingle,mu)
		Bpar[0,j,i] = B[i] * np.cos(thetas[j]*np.pi/180.0)
		Bperp[j,i] = B[i] * np.sin(thetas[j]*np.pi/180.0)
		SV[0,j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)
		SQ[j,i] = -scipy.integrate.simps(stokes[1,:] * pQ, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pQ, x=wavelength)
		
		modelSingle = np.asarray([B[i], 180.0-thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
		stokes = s.synth(modelSingle,mu)
		Bpar[1,j,i] = B[i] * np.cos((180.0-thetas[j])*np.pi/180.0)
		SV[1,j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)		
pl.close('all')
f, ax = pl.subplots()
for j in range(4):
	ax.plot(Bpar[0,j,:], SV[0,j,:], color=pal[j], label=r'$\theta_B=${0}'.format(thetas[j]))
	ax.plot(Bpar[1,j,:], SV[1,j,:], color=pal[j])
	
ax.set_xlabel('B$_\parallel$')
ax.set_ylabel('S$_V$')
ax.set_xlim([-3500,3500])
pl.legend()
f, ax = pl.subplots()
for j in range(4):
	ax.plot(Bperp[j,:], SQ[j,:], color=pal[j], label=r'$\theta_B=${0}'.format(thetas[j]))	
	
ax.set_xlabel('B$_\perp$')
ax.set_ylabel('S$_Q$')
ax.set_xlim([0,3500])
pl.legend(loc='upper left')
Explanation: Longitudinal magnetograph
<a id='longMagnetograph'></a>
End of explanation
lambda0 = 5250.21
JUp = 0.0
JLow = 1.0
gUp = 0.0
gLow = 3.0
lambdaStart = 5249.5
lambdaStep = 0.015
nLambda = 100
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep
lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])
s = pymilne.milne(nLambda, lineInfo)
BField = 100.0
BTheta = 0.0
BChi = 0.0
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 12.0
mu = 1.0
VDop = 0.032
kl = 6.1
nB = 100
B = np.linspace(0.0,3500.0,nB)
Bpar = np.zeros((4,nB))
BEstimated = np.zeros((4,nB))
pal = sn.color_palette()
thetas = [15, 30, 60, 75]
for j in range(4):
	for i in range(nB):
		modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
		stokes = s.synth(modelSingle,mu)
		
		Ic = stokes[0,0]
		
		Iplus = 0.5*(stokes[0,:] + stokes[3,:])
		Iminus = 0.5*(stokes[0,:] - stokes[3,:])
		
		lambdaPlus = scipy.integrate.simps((0.5*Ic - Iplus) * wavelength, x=wavelength) / scipy.integrate.simps((0.5*Ic - Iplus), x=wavelength)
		lambdaMinus = scipy.integrate.simps((0.5*Ic - Iminus) * wavelength, x=wavelength) / scipy.integrate.simps((0.5*Ic - Iminus), x=wavelength)
		
		BEstimated[j,i] = 1.071e9 / (3.0 * 5250.21**2) * (lambdaPlus-lambdaMinus)*1e3
		Bpar[j,i] = B[i] * np.cos(thetas[j]*np.pi/180.0)		
pl.close('all')
f, ax = pl.subplots()
for j in range(4):
	ax.plot(Bpar[j,:], (BEstimated[j,:]-Bpar[j,:]) / Bpar[j,:], color=pal[j], label=r'$\theta_B=${0}'.format(thetas[j]))	
	
ax.set_xlabel('B$_\parallel$')
ax.set_ylabel('Relative error')
ax.set_xlim([0,3500])
pl.legend()
Explanation: Center-of-gravity
<a id='cog'></a>
End of explanation
thetaReal = np.linspace(0.0,180.0,200)
x = np.asarray([1.0,0.1,0.01])
f, ax = pl.subplots()
tanTheta = np.tan(thetaReal*np.pi/180.0)
ind = np.where(tanTheta < 0.0)
for i in range(3):
    out = np.arctan2(tanTheta, np.sqrt(x[i]))
    out[ind] += np.pi
    ax.plot(thetaReal, out*180.0/np.pi, label='x={0}'.format(x[i]))
ax.set_xlabel(r'$\theta_m$')
ax.set_ylabel(r'$\theta_\mathrm{app}$')
pl.legend(loc='upper left')
Explanation: Influence of filling factor on inclination of magnetic field
<a id='inclinationFilling'></a>
End of explanation
# First line
lambda0 = 6302.5
JUp = 1.0
JLow = 1.0
gUp = 2.5
gLow = 0.0
lambdaStart = 6301.8
lambdaStep = 0.03
nLambda = 50
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep
lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])
s = pymilne.milne(nLambda, lineInfo)
VMac = 0.0
damping = 0.04
B0 = 1.0
B1 = 5.0
mu = 1.0
VDop = 0.040
kl = 6.75
pV = np.exp(-(wavelength-lambda0-60e-3)**2 / 75e-3**2)
nB = 100
B = np.linspace(0.0,3500.0,nB)
SV1 = np.zeros((4,nB))
SV2 = np.zeros((4,nB))
pal = sn.color_palette()
thetas = [0, 30, 60, 85]
for j in range(4):
    for i in range(nB):
        modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
        stokes = s.synth(modelSingle,mu)
        SV1[j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)
        
# Second line
lambda0 = 6301.5
JUp = 1.0
JLow = 1.0
gUp = 1.5
gLow = 0.0
lambdaStart = 6300.8
lambdaStep = 0.03
nLambda = 50
wavelength = lambdaStart + np.arange(nLambda) * lambdaStep
lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])
s = pymilne.milne(nLambda, lineInfo)
pV = np.exp(-(wavelength-lambda0-60e-3)**2 / 75e-3**2)
thetas = [0, 30, 60, 85]
for j in range(4):
    for i in range(nB):
        modelSingle = np.asarray([B[i], thetas[j], 0.0, VMac, damping, B0, B1, VDop, kl])
        stokes = s.synth(modelSingle,mu)
        SV2[j,i] = scipy.integrate.simps(stokes[3,:] * pV, x=wavelength) / scipy.integrate.simps(stokes[0,:] * pV, x=wavelength)
        
f, ax = pl.subplots()
for i in range(4):
    ax.plot(B, SV1[i,:] / SV2[i,:])
for i in range(4):
    ax.plot(B, SV1[i,:] / SV2[i,:], color=pal[i], label=r'$\theta_B$={0}'.format(thetas[i]))
ax.set_xlabel('B [G]')
ax.set_ylabel('Line ratio')
ax.axhline(2.5/1.5, color=pal[4])
pl.legend(loc='lower left')
Explanation: Line ratio technique
<a id='lineRatio'></a>
End of explanation
<END_TASK>