Skip to content

Commit 37e0d76

Browse files
RevathiJambunathanRevathiJambunathan
committed
add oblique angle Chi based on Eq 2.7, 2.8, 2.9 in Jeromi Petris 2016 article
1 parent 2bd3a5b commit 37e0d76

File tree

3 files changed

+177
-82
lines changed

3 files changed

+177
-82
lines changed

‎Source/Particles/PhysicalParticleContainer.cpp

+6-4
Original file line numberDiff line numberDiff line change
@@ -2320,6 +2320,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
23202320
amrex::Real Bstar_data = Pulsar::m_B_star;
23212321
amrex::Real Rstar_data = Pulsar::m_R_star;
23222322
amrex::Real dRstar_data = Pulsar::m_dR_star;
2323+
amrex::Real chi_data = Pulsar::m_Chi;
23232324
amrex::Real corotatingE_maxradius_data = Pulsar::m_corotatingE_maxradius;
23242325
int E_external_monopole_data = Pulsar::m_do_E_external_monopole;
23252326
int AddExternalMonopoleOnly = Pulsar::m_AddExternalMonopoleOnly;
@@ -2425,7 +2426,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
24252426
amrex::ParticleReal r_p, theta_p, phi_p;
24262427
Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr,
24272428
r_p, theta_p, phi_p);
2428-
Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p,
2429+
Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, chi_data,
24292430
AddExternalMonopoleOnly,
24302431
AddPulsarVacuumEFields,
24312432
AddBDipoleOutsideRstar,
@@ -2436,7 +2437,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
24362437
Bstar_data, Rstar_data, dRstar_data,
24372438
Exp, Eyp, Ezp, Bxp, Byp, Bzp);
24382439
if (use_theoreticalEB == 1) {
2439-
Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p,
2440+
Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p, chi_data,
24402441
theory_max_rstar,
24412442
corotatingE_maxradius_data,
24422443
E_external_monopole_data,
@@ -2747,6 +2748,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
27472748
amrex::Real Bstar_data = Pulsar::m_B_star;
27482749
amrex::Real Rstar_data = Pulsar::m_R_star;
27492750
amrex::Real dRstar_data = Pulsar::m_dR_star;
2751+
amrex::Real chi_data = Pulsar::m_Chi;
27502752
amrex::Real corotatingE_maxradius_data = Pulsar::m_corotatingE_maxradius;
27512753
int E_external_monopole_data = Pulsar::m_do_E_external_monopole;
27522754
int AddExternalMonopoleOnly = Pulsar::m_AddExternalMonopoleOnly;
@@ -2878,7 +2880,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
28782880
amrex::ParticleReal r_p, theta_p, phi_p;
28792881
Pulsar::ConvertCartesianToSphericalCoord( xp, yp, zp, center_star_arr,
28802882
r_p, theta_p, phi_p);
2881-
Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p,
2883+
Pulsar::getExternalEBOnParticle(r_p, theta_p, phi_p, chi_data,
28822884
AddExternalMonopoleOnly,
28832885
AddPulsarVacuumEFields,
28842886
AddBDipoleOutsideRstar,
@@ -2889,7 +2891,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
28892891
Bstar_data, Rstar_data, dRstar_data,
28902892
Exp, Eyp, Ezp, Bxp, Byp, Bzp);
28912893
if (use_theoreticalEB == 1) {
2892-
Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p,
2894+
Pulsar::EnforceTheoreticalEBOnParticle(r_p, theta_p, phi_p, chi_data,
28932895
theory_max_rstar,
28942896
corotatingE_maxradius_data,
28952897
E_external_monopole_data,

‎Source/Particles/PulsarParameters.H

+97-39
Original file line numberDiff line numberDiff line change
@@ -58,34 +58,49 @@ public:
5858

5959
AMREX_GPU_HOST_DEVICE AMREX_INLINE
6060
static void CorotatingEfieldSpherical (amrex::Real const r, amrex::Real const theta,
61-
amrex::Real const phi, amrex::Real const time,
61+
amrex::Real const phi, amrex::Real const chi,
62+
amrex::Real const time,
6263
amrex::Real const omega_star_data,
6364
amrex::Real const ramp_omega_time_data,
6465
amrex::Real const Bstar,
6566
amrex::Real const Rstar,
6667
amrex::Real const dRstar,
6768
amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi)
6869
{
69-
amrex::ignore_unused(phi);
7070
amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data);
71+
// Polar angle
7172
amrex::Real c_theta = std::cos(theta);
7273
amrex::Real s_theta = std::sin(theta);
74+
// Oblique angle
75+
amrex::Real c_chi = std::cos(chi);
76+
amrex::Real s_chi = std::sin(chi);
77+
// Instantaneous phase, psi = phi - Omega*t (phi is the azimuthal angle)
78+
// Refer to Pg. 6 of Jeromi Petri's 2016 paper
79+
amrex::Real psi = phi - omega*time;
80+
amrex::Real c_psi = std::cos(psi);
81+
// Ratio : Rstar/r
7382
amrex::Real r_ratio;
7483
if (r > 0) {
7584
r_ratio = Rstar/r;
7685
} else {
7786
r_ratio = Rstar/(dRstar*0.5);
7887
}
7988
amrex::Real r2 = r_ratio * r_ratio;
80-
// Michel and Li -- eq 14 , 15
81-
Er = Bstar * omega * r2 * Rstar * s_theta * s_theta;
82-
Etheta = -Bstar * omega * r2 * Rstar * 2.0 * s_theta * c_theta;
83-
Ephi = 0.0; // aligned magnetic and rotation axis
89+
// Corotating electric field inside the pulsar that corresponds to dipole magnetic field
90+
// Eq 2.9 (a-c) Jeromi Petri 2016
91+
Er = Bstar * omega * r2 * Rstar *
92+
( c_chi * s_theta * s_theta
93+
- s_chi * c_theta * s_theta * c_psi);
94+
Etheta = -Bstar * omega * r2 * Rstar *
95+
( c_chi * 2._rt * s_theta * c_theta
96+
+ s_chi * 2._rt * s_theta * s_theta * c_psi);
97+
Ephi = 0._rt;
8498
}
8599

86100
AMREX_GPU_HOST_DEVICE AMREX_INLINE
87101
static void ExternalEFieldSpherical (amrex::Real const r, amrex::Real const theta,
88-
amrex::Real const phi, amrex::Real const time,
102+
amrex::Real const phi, amrex::Real const chi,
103+
amrex::Real const time,
89104
amrex::Real const omega_star_data,
90105
amrex::Real const ramp_omega_time_data,
91106
amrex::Real const Bstar, amrex::Real const Rstar,
@@ -94,38 +109,59 @@ public:
94109
int const ApplyCorotatingEField,
95110
amrex::Real &Er, amrex::Real &Etheta, amrex::Real &Ephi)
96111
{
97-
amrex::ignore_unused(phi, corotatingE_maxradius);
112+
amrex::ignore_unused(corotatingE_maxradius);
98113
amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data);
114+
// Polar angle
99115
amrex::Real c_theta = std::cos(theta);
100116
amrex::Real s_theta = std::sin(theta);
117+
// Oblique angle
118+
amrex::Real c_chi = std::cos(chi);
119+
amrex::Real s_chi = std::sin(chi);
120+
// Instantaneous phase, psi = phi - Omega*t (phi is the azimuthal angle)
121+
// Refer to Pg. 6 of Jeromi Petri's 2016 paper
122+
amrex::Real psi = phi - omega*time;
123+
amrex::Real c_psi = std::cos(psi);
124+
amrex::Real s_psi = std::sin(psi);
125+
101126
amrex::Real r_ratio = Rstar/r;
102127
// inside pulsar
103-
//if (r <= corotatingE_maxradius) {
104128
if (ApplyCorotatingEField == 1) {
105129
amrex::Real r2 = r_ratio * r_ratio;
106-
// Michel and Li -- eq 14 , 15
107-
Er = Bstar * omega * r2 * Rstar * s_theta * s_theta;
108-
Etheta = -Bstar * omega * r2 * Rstar * 2.0 * s_theta * c_theta;
109-
Ephi = 0.0; // aligned magnetic and rotation axis
130+
// Corotating electric field inside the pulsar that corresponds to dipole magnetic field
131+
// Eq 2.9 (a-c) Jeromi Petri 2016
132+
Er = Bstar * omega * r2 * Rstar *
133+
( c_chi * s_theta * s_theta
134+
- s_chi * c_theta * s_theta * c_psi);
135+
Etheta = -Bstar * omega * r2 * Rstar *
136+
( c_chi * 2._rt * s_theta * c_theta
137+
+ s_chi * 2._rt * s_theta * s_theta * c_psi);
138+
Ephi = 0._rt;
110139
}
111140

112141
// outside pulsar
113-
//if (r > corotatingE_maxradius) {
114142
if (ApplyCorotatingEField == 0) {
115-
amrex::Real r4 = r_ratio*r_ratio*r_ratio*r_ratio;
116-
// Taking derivative of phi given in eq 30 of Michel and Li
117-
Er = Bstar * omega * Rstar * r4 * (1.0-3.0*c_theta*c_theta);
143+
amrex::Real r2 = r_ratio * r_ratio;
144+
amrex::Real r4 = r2 * r2;
145+
// Equations 2.8(a) - 2.8(c)
146+
Er = Bstar * omega * Rstar * r4 *
147+
( c_chi * (1._rt - 3._rt * c_theta * c_theta)
148+
- 3._rt * s_chi * c_theta * s_theta * c_psi );
118149
if (Eexternal_monopole == 1) {
119-
Er += (2.0/3.0) * omega * Bstar * Rstar * r_ratio * r_ratio;
150+
// Monopole term from central charge = Qc/(4*pi*eps0*r^2)
151+
// For dipolar magnetic field, Qc = 8pi/3 * eps0*Omega*Bstar*Rstar^3*c_chi
152+
Er += (2._rt/3._rt) * omega * Bstar * Rstar * r2 * c_chi;
120153
}
121-
Etheta = (-1.0) * Bstar * omega * Rstar * r4 * (2.0*s_theta*c_theta);
122-
Ephi = 0.0;
154+
Etheta = Bstar * omega * Rstar *
155+
( r2 * s_chi * (r2*std::cos(2._rt*theta) - 1._rt) * c_psi
156+
- r4 * c_chi * std::sin(2._rt*theta) );
157+
Ephi = Bstar * omega * Rstar * r2 * (1._rt - r2) * s_chi * c_theta * s_psi;
123158
}
124159
}
125160

126161
AMREX_GPU_HOST_DEVICE AMREX_INLINE
127162
static void ExternalEMonopoleSpherical (amrex::Real const r, amrex::Real const theta,
128-
amrex::Real const phi, amrex::Real const time,
163+
amrex::Real const phi, amrex::Real const chi,
164+
amrex::Real const time,
129165
amrex::Real const omega_star_data,
130166
amrex::Real const ramp_omega_time_data,
131167
amrex::Real const Bstar, amrex::Real const Rstar,
@@ -134,39 +170,55 @@ public:
134170
amrex::ignore_unused(phi, theta);
135171
amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data);
136172
amrex::Real r_ratio = Rstar/r;
173+
amrex::Real r2 = r_ratio * r_ratio;
174+
amrex::Real c_chi = std::cos(chi);
137175

138-
Er = (2.0/3.0) * omega * Bstar * Rstar * r_ratio * r_ratio;
176+
Er = (2._rt/3._rt) * omega * Bstar * Rstar * r2 * c_chi;
139177
Etheta = 0.;
140178
Ephi = 0.;
141179

142180
}
143181

144182
AMREX_GPU_HOST_DEVICE AMREX_INLINE
145183
static void ExternalBFieldSpherical (amrex::Real const r, amrex::Real const theta,
146-
amrex::Real const phi, amrex::Real const time,
184+
amrex::Real const phi, amrex::Real const chi,
185+
amrex::Real const time,
186+
amrex::Real const omega_star_data,
187+
amrex::Real const ramp_omega_time_data,
147188
amrex::Real const Bstar, amrex::Real const Rstar,
148189
amrex::Real const dRstar,
149190
amrex::Real &Br, amrex::Real &Btheta, amrex::Real &Bphi)
150191
{
151-
amrex::ignore_unused(phi, time);
192+
// Polar angle
152193
amrex::Real c_theta = std::cos(theta);
153194
amrex::Real s_theta = std::sin(theta);
195+
// Oblique angle
196+
amrex::Real c_chi = std::cos(chi);
197+
amrex::Real s_chi = std::sin(chi);
198+
// Instantaneous phase, psi = phi - Omega*t (phi is the azimuthal angle)
199+
// Refer to Pg. 6 of Jeromi Petri's 2016 paper
200+
amrex::Real omega = Omega(omega_star_data, time, ramp_omega_time_data);
201+
amrex::Real psi = phi - omega*time;
202+
amrex::Real c_psi = std::cos(psi);
203+
amrex::Real s_psi = std::sin(psi);
204+
154205
amrex::Real r_ratio;
155206
if (r > 0) {
156207
r_ratio = Rstar/r;
157208
} else {
158209
r_ratio = Rstar/(dRstar*0.5);
159210
}
160211
amrex::Real r3 = r_ratio*r_ratio*r_ratio;
161-
// Michel and Li -- eq 14 and 15 from michel and li
162-
// dipole B field inside and outside the pulsar
163-
Br = 2.0*Bstar*r3*c_theta;
164-
Btheta = Bstar*r3*s_theta;
165-
Bphi = 0.0;
212+
// The full dipole magnetic field as f(theta, phi, omega, t, and Chi)
213+
// Eqs 2.7(a) - 2.7(c) from Jeromi Petri 2016 paper
214+
Br = 2._rt * Bstar * r3 * ( c_chi * c_theta + s_chi * s_theta * c_psi);
215+
Btheta = Bstar * r3 * ( c_chi * s_theta - s_chi * c_theta * c_psi);
216+
Bphi = Bstar * r3 * s_chi * s_psi;
166217
}
167218

168219
AMREX_GPU_HOST_DEVICE AMREX_INLINE
169-
static void getExternalEBOnParticle( amrex::Real const r, amrex::Real const theta, amrex::Real const phi,
220+
static void getExternalEBOnParticle( amrex::Real const r, amrex::Real const theta,
221+
amrex::Real const phi, amrex::Real const chi,
170222
const int AddExternalMonopoleOnly, const int AddPulsarVacuumEFields,
171223
const int AddBDipoleOutsideRstar, const int AddPulsarVacuumBFields,
172224
amrex::Real const corotatingE_maxradius,
@@ -184,7 +236,7 @@ public:
184236
amrex::Real Er = 0.0_rt;
185237
amrex::Real Etheta = 0.0_rt;
186238
amrex::Real Ephi = 0.0_rt;
187-
ExternalEMonopoleSpherical(r, theta, phi, cur_time, omega_star,
239+
ExternalEMonopoleSpherical(r, theta, phi, chi, cur_time, omega_star,
188240
ramp_omega_time, Bstar, Rstar,
189241
Er, Etheta, Ephi);
190242
amrex::Real Ex_monopole = 0._rt;
@@ -206,7 +258,7 @@ public:
206258
amrex::Real Ephi = 0.0_rt;
207259
int ApplyCorotatingEField = 0;
208260
if (r <= corotatingE_maxradius) ApplyCorotatingEField = 1;
209-
ExternalEFieldSpherical(r, theta, phi, cur_time, omega_star,
261+
ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star,
210262
ramp_omega_time, Bstar, Rstar,
211263
corotatingE_maxradius,
212264
E_external_monopole,
@@ -229,8 +281,10 @@ public:
229281
amrex::Real Br = 0._rt;
230282
amrex::Real Btheta = 0._rt;
231283
amrex::Real Bphi = 0._rt;
232-
ExternalBFieldSpherical (r, theta, phi, cur_time, Bstar, Rstar, dR_star,
233-
Br, Btheta, Bphi);
284+
ExternalBFieldSpherical(r, theta, phi, chi, cur_time,
285+
omega_star, ramp_omega_time,
286+
Bstar, Rstar, dR_star,
287+
Br, Btheta, Bphi);
234288
amrex::Real Bx_dipole = 0._rt;
235289
amrex::Real By_dipole = 0._rt;
236290
amrex::Real Bz_dipole = 0._rt;
@@ -246,8 +300,10 @@ public:
246300
amrex::Real Br = 0._rt;
247301
amrex::Real Btheta = 0._rt;
248302
amrex::Real Bphi = 0._rt;
249-
ExternalBFieldSpherical (r, theta, phi, cur_time, Bstar, Rstar, dR_star,
250-
Br, Btheta, Bphi);
303+
ExternalBFieldSpherical(r, theta, phi, chi, cur_time,
304+
omega_star, ramp_omega_time,
305+
Bstar, Rstar, dR_star,
306+
Br, Btheta, Bphi);
251307
amrex::Real Bx_dipole = 0._rt;
252308
amrex::Real By_dipole = 0._rt;
253309
amrex::Real Bz_dipole = 0._rt;
@@ -262,7 +318,8 @@ public:
262318

263319

264320
AMREX_GPU_HOST_DEVICE AMREX_INLINE
265-
static void EnforceTheoreticalEBOnParticle( amrex::Real const r, amrex::Real const theta, amrex::Real const phi,
321+
static void EnforceTheoreticalEBOnParticle( amrex::Real const r, amrex::Real const theta,
322+
amrex::Real const phi, amrex::Real const chi,
266323
amrex::Real const theory_max_radius,
267324
amrex::Real const corotatingE_maxradius,
268325
const int E_external_monopole,
@@ -285,7 +342,7 @@ public:
285342
amrex::Real Ex_theory = 0._rt;
286343
amrex::Real Ey_theory = 0._rt;
287344
amrex::Real Ez_theory = 0._rt;
288-
ExternalEFieldSpherical(r, theta, phi, cur_time, omega_star, ramp_omega_time,
345+
ExternalEFieldSpherical(r, theta, phi, chi, cur_time, omega_star, ramp_omega_time,
289346
Bstar, Rstar, corotatingE_maxradius, E_external_monopole,
290347
ApplyCorotatingEField, Er, Etheta, Ephi);
291348
// Compute theoretical Bfield in spherical coordinates
@@ -295,7 +352,8 @@ public:
295352
amrex::Real Bx_theory = 0._rt;
296353
amrex::Real By_theory = 0._rt;
297354
amrex::Real Bz_theory = 0._rt;
298-
ExternalBFieldSpherical(r, theta, phi, cur_time, Bstar, Rstar, dR_star, Br, Btheta, Bphi);
355+
ExternalBFieldSpherical(r, theta, phi, chi, cur_time, omega_star, ramp_omega_time,
356+
Bstar, Rstar, dR_star, Br, Btheta, Bphi);
299357
// Convert Efield from spherical to Cartesian
300358
ConvertSphericalToCartesianXComponent(Er, Etheta, Ephi, r, theta, phi, Ex_theory);
301359
ConvertSphericalToCartesianYComponent(Er, Etheta, Ephi, r, theta, phi, Ey_theory);

0 commit comments

Comments
 (0)