Notes

💡  These are the functions described and used in van Leeuwen & Timokhin 2012 (LT12). 
     You are encouraged to use them for your own research; in that case kindly cite the paper.

💡  Sections in purple will be evaluated automatically, and contain the functions to calculate the potential drop
💡  Only sections with black titles needs to be evaluated: insert your own data there; the evaluation order is important

Setup

SetDirectory[NotebookDirectory[]]

Pulsar Parameters

ξ1PSR=1;
ΔξPSR=0.2;

Solution for coefficients in the expression for V(ξ) - auxilary functions

Coefficients a,b are determined from dV/dξ = 2 aξ + b = ξ((Subscript[ω, D])/180°) for bands 1,2;
 v1, v2  are  ((Subscript[ω, D 1,2])/180°)

ab[v1_,v2_,Δξ_,ξ1_]:=
Module[{x1,x2},
Solve[{
v1 x1 == 2a x1 +b,
v2 x2 == 2a x2 +b
}/.{x2->ξ1-Δξ,x1->ξ1},
{a,b}]]

Expression for the electric potential (assuming V((Subscript[ξ, 1]))=0)

V[ξ_,v1_,v2_,Δξ_,ξ1_]:=a ((ξ^2)-((ξ1)^2))+b (ξ-ξ1) /.ab[v1,v2,Δξ,ξ1][[1]]

Expression for the electric potential for PSR B0826-34 
vdeg1, vdeg2 are  (Subscript[ω, D 1,2]) - drift speeds at points 1, 2 in [degrees/Period]

VPSR0826[ξ_,vdeg1_,vdeg2_]:=Module[{v1=(vdeg1/180),v2=(vdeg2/180)},V[ξ,v1,v2,ΔξPSR,ξ1PSR]]

Check for consistency

D[V[ξ,v1,v2,Δξ,ξ1],{ξ}]/.ξ->ξ1//Simplify

D[V[ξ,v1,v2,Δξ,ξ1],{ξ}]/.ξ->ξ1-Δξ//Simplify

V[ξ,v1,v2,Δξ,1]//Simplify

Plotting function for the electric potential: vdeg1, vdeg2 -- drift speeds at points 1,2 in [degrees/Period]

This function plots potential V(ξ) for ξ=[ξ1PSR-ΔξPSR, ξ1PSR0] if you give drift speeds (Subscript[v, D 1]),(Subscript[v, D 2])in [ degrees/Period ]

plotVPSR0826[vdeg1_,vdeg2_]:=
Module[ {v1=(vdeg1/180),v2=(vdeg2/180)} ,
Plot[ Evaluate[ V[ξ,v1,v2,ΔξPSR,ξ1PSR] ],{ξ,ξ1PSR-ΔξPSR,ξ1PSR} ]]

Example:

plotVPSR0826[2,-2]

plotVPSR0826[2,2]

plotVPSR0826[2,-4]

Maximum Potential drop

Local extremum position

ξex[v1_,v2_,Δξ_,ξ1_]:= -(b/(2 a))/.ab[v1,v2,Δξ,ξ1][[1]]

general expression for the maximum potential drop

ΔVmax[v1_,v2_,Δξ_,ξ1_]:= 
Module[
{ξ2=ξ1-Δξ},
If[ ξex[v1,v2,Δξ,ξ1]>ξ1 || ξex[v1,v2,Δξ,ξ1]<ξ2,

V[ξ2,v1,v2,Δξ,ξ1],

Module[
{Ve1=-(a ((ξ1)^2) + b ξ1 +((b^2)/(4 a)))/.ab[v1,v2,Δξ,ξ1][[1]],
Ve2=-(a ((ξ2)^2) + b ξ2 +((b^2)/(4 a)))/.ab[v1,v2,Δξ,ξ1][[1]]},
 If[Abs[Ve1]≥ Abs[Ve2], Ve1, Ve2]]]]

Maximum Potential Drop for PSR B 0826-34

ΔVmaxPSR0826[vdeg1_,vdeg2_]:=Module[{v1=(vdeg1/180),v2=(vdeg2/180)},ΔVmax[v1,v2,ΔξPSR,ξ1PSR]]

ΔVmaxPSR0826[2,-2]

List of pulses for plotting Potential drop, LT12 paper.

selectListVExample={2,12,21, 24, 30};

Drift rates of PSR B 0826-34 as used in LT12 paper. Insert own data here.

Prepare data

Uncomment the line below to read your own data from a text file

(* driftRatesList=ReadList["YourDataHere.txt",Number,RecordLists->True] *)

Or use the data from the LT12 paper:

driftRatesList={{40.`,0.630154`,-0.674995`},{50.`,0.587671`,-0.816604`},{60.`,0.665555`,-0.495626`},{70.`,0.821324`,0.231292`},{80.`,0.38942`,0.665555`},{90.`,0.261974`,0.56643`},{100.`,-0.261974`,0.066083`},{110.`,0.247813`,0.0330421`},{120.`,0.0849643`,-0.0849645`},{160.`,1.16118`,1.416076`},{170.`,0.644315`,0.165208`},{180.`,0.247813`,-0.174648`},{190.`,-0.453142`,-0.0708044`},{200.`,0.74344`,0.0708038`},{250.`,1.16`,1.033736`},{280.`,0.9`,0.2926552`},{320.`,0.6`,-0.46`},{330.`,-0.509787`,-0.25`},{380.`,0.162848`,0.127448`},{390.`,-0.0424827`,-0.297377`},{400.`,0.0637275`,-0.771761`},{410.`,-0.686799`,-0.58059`},{590.`,0.5522702`,1.274467`},{620.`,1.26`,-0.1`},{630.`,-0.1557663`,-0.7788323`},{700.`,-0.332778`,0.302099`},{710.`,-0.0566422`,-0.354021`},{720.`,-0.113287`,-0.292655`},{730.`,-0.375262`,-0.438984`},{740.`,-0.502703`,-0.462583`},{750.`,-0.276135`,-0.231294`},{760.`,-0.304453`,-0.542828`},{770.`,-0.1274498`,-0.745798`},{780.`,0.240731`,0.44606474`},{790.`,-0.070803069`,0.26197146`}}

driftRatesList//TableForm

Plotting the drift rates of PSR B 0826-34 (Fig. 7 in LT12)

Store specific portions of data in variables

driftRatesList: (Subscript[ξ, 2]) <==> Column 2  : Left band, close to the center, in our notations at (Subscript[ξ, 2])
driftRatesList: (Subscript[ξ, 1]) <==> Column 3  : Right band, close to the boundary, in our notations at (Subscript[ξ, 1])

pulseN=driftRatesList[[All,1]]

pulseDrift1=driftRatesList[[All,{1,3}]]

pulseDrift2=driftRatesList[[All,{1,2}]]

driftRatesList[[selectListVExample,1]]

Make Plot for drift speed

driftPlotBase[vmin_,vmax_]:=
ListLinePlot[
{
pulseDrift1[[1;;9]],pulseDrift2[[1;;9]],
pulseDrift1[[10;;14]],pulseDrift2[[10;;14]],
pulseDrift1[[15;;15]],pulseDrift2[[15;;15]],
pulseDrift1[[16;;16]],pulseDrift2[[16;;16]],
pulseDrift1[[17;;18]],pulseDrift2[[17;;18]],
pulseDrift1[[19;;22]],pulseDrift2[[19;;22]],
pulseDrift1[[23;;23]],pulseDrift2[[23;;23]],
pulseDrift1[[24;;25]],pulseDrift2[[24;;25]],
pulseDrift1[[26;;35]],pulseDrift2[[26;;35]]
},
PlotRange->{{0,All},{vmin,vmax}},
PlotMarkers->{Automatic,5} ,
PlotStyle->
{{AbsoluteThickness[0.25],Black},
{AbsoluteDashing[3],AbsoluteThickness[0.25],Black}},
Frame->True,
FrameLabel->
{{Style["(Subscript[ω, D]) [°/Period]",FontFamily->"Times",FontSlant->Plain,FontSize->8],None},{Style["pulse number",FontFamily->"Times",FontSlant->Plain,FontSize->8] ,None}},
FrameTicksStyle->Directive[7]
]

driftPlot[selectList_]:=
Module[{xmin=0, ymin=-1.3,ymax=1.9, yselectlabel=-1.,selectlabel={a,b,c,d,e,f,g,h}},
Show[
{driftPlotBase[ymin,ymax],
Graphics[
{AbsoluteDashing[{0.05,2}],
AbsoluteThickness[0.25],
Table[Line[ {{pulseN[[i]],ymin},{pulseN[[i]],ymax}} ],{i,selectList}]}],
Graphics[
{FontSize->8,
Table[Text[selectlabel[[i]], {pulseN[[selectList[[i]]]]+20,yselectlabel}],{i,Length[selectList]}]}]} ,
ImageSize->72*3.5]]

paperDriftPlot=driftPlot[selectListVExample]

Export["drift_speed.pdf",paperDriftPlot]

Maximum potential drop variation (Fig. 8 in LT12)

Prepare data

driftRatesList: (Subscript[ξ, 2]) <==> Column 2  : Left band, close to the center, in our notations at (Subscript[ξ, 2])
driftRatesList: (Subscript[ξ, 1]) <==> Column 3  : Right band, close to the boundary, in our notations at (Subscript[ξ, 1])

driftRates12=driftRatesList[[All,{3,2}]]

dataΔVmax=ΔVmaxPSR0826[#[[1]],#[[2]]]&/@driftRates12

pulseΔVmax=Thread[({#1,#2}&)[driftRatesList[[All,1]],dataΔVmax]]

Make plot for maximum potential drop variation

VmaxPlotBase[dVmin_,dVmax_]:=
ListLinePlot[
{
pulseΔVmax[[1;;9]],
pulseΔVmax[[10;;14]],
pulseΔVmax[[15;;15]],
pulseΔVmax[[16;;16]],
pulseΔVmax[[17;;18]],
pulseΔVmax[[19;;22]],
pulseΔVmax[[23;;23]],
pulseΔVmax[[24;;25]],
pulseΔVmax[[26;;35]]
},
PlotRange->{{0,All},{dVmin,dVmax}},
PlotMarkers->{Automatic,5} ,
PlotStyle->{{AbsoluteThickness[0.25],Black} ,{AbsoluteThickness[0.25],Black} },
Frame->True,
FrameLabel->
{{Style["(Subscript[ΔV, max])",FontFamily->"Times",FontSlant->Plain,FontSize->8],None},{Style["pulse number",FontFamily->"Times",FontSlant->Plain,FontSize->8] ,None}},
FrameTicksStyle->Directive[7]
]

VmaxPlot[selectList_]:=Module[{ymin=-0.0015,ymax=0.0011, yselectlabel=0.00095,selectlabel={a,b,c,d,e,f,g,h}},
Show[
{VmaxPlotBase[ymin,ymax],
Graphics[
{AbsoluteDashing[{0.05,2}],
AbsoluteThickness[0.25],
Table[Line[ {{pulseN[[i]],ymin},{pulseN[[i]],ymax}} ],{i,selectList}]}],
Graphics[
{FontSize->8,
Table[Text[selectlabel[[i]], {pulseN[[selectList[[i]]]]+20,yselectlabel}],{i,Length[selectList]}]}]} ,
ImageSize->72*3.5]]

paperVmaxPlot=VmaxPlot[selectListVExample]

Export["vmax.pdf",paperVmaxPlot]

Example of Potential drop variation (Fig. 9 in LT12)

Prepare data

driftRatesList: (Subscript[ξ, 2]) <==> Column 2  : Left band, close to the center, in our notations at (Subscript[ξ, 2])
driftRatesList: (Subscript[ξ, 1]) <==> Column 3  : Right band, close to the boundary, in our notations at (Subscript[ξ, 1])

drift1=driftRatesList[[All,3]]

drift2=driftRatesList[[All,2]]

Make Plot with example of potential drop variations

VExamplePlotBase[selectList_,{{xmin_,xmax_},{Vmin_,Vmax_}}]:=
Plot[
Evaluate[ Table[ VPSR0826[ξ,drift1[[i]],drift2[[i]]],{i,selectList}] ],
{ξ,ξ1PSR-ΔξPSR,ξ1PSR},
PlotRange->{{xmin,xmax},{Vmin,Vmax}},
PlotStyle->{{AbsoluteThickness[0.5],Black}} ,
Axes->False,
Frame->True,
FrameLabel->
{{Style["V(ξ)",FontFamily->"Times",FontSlant->Plain,FontSize->8],None},{Style["ξ",FontFamily->"Times",FontSlant->Plain,FontSize->8] ,None}},
FrameTicksStyle->Directive[7]]

VExamplePlot[selectList_]:=
Module[
{xmin=0.78,xmax=1.005,ymin=-0.0008,ymax=0.0008,
xselectlabel=0.79,
selectlabel={a,b,c,d,e,f,g,h}},
Show[
{VExamplePlotBase[selectList,{{xmin,xmax},{ymin,ymax}}],
Graphics[
{FontSize->8,
Table[
Text[selectlabel[[i]], 
{xselectlabel, VPSR0826[ξ1PSR-ΔξPSR,drift1[[ selectList[[i]] ]],drift2[[ selectList[[i]]] ]] }],
{i,Length[selectList]}]}] 
},
ImageSize->72*3.5]]

paperVExamplePlot=VExamplePlot[ selectListVExample ]

Export["v_example.pdf",paperVExamplePlot ]