Repository: fatiando/harmonica
Pull request #414 (2023-05-19 21:38:38)

Magnetic field of dipoles in Cartesian coordinates
Add forward modelling functions for the magnetic field of dipoles in Cartesian coordinates: `dipole_magnetic` and `dipole_magnetic_component`. The first one computes the three components of the magnetic field, while the latter computes only one of them. The forward modelling is done using the Choclo functions for dipoles.

**Relevant issues/PRs:**

Part of #406




Patch:
diff --git a/doc/_static/figures/dipoles.svg b/doc/_static/figures/dipoles.svg
new file mode 100644
index 0000000..f3d0b3f
--- /dev/null
+++ b/doc/_static/figures/dipoles.svg
@@ -0,0 +1,625 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!-- Created with Inkscape (http://www.inkscape.org/) -->
+
+<svg
+   width="100"
+   height="100"
+   viewBox="0 0 26.458334 26.458334"
+   version="1.1"
+   id="svg344"
+   inkscape:version="1.2.2 (b0a8486541, 2022-12-01)"
+   sodipodi:docname="dipoles.svg"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:svg="http://www.w3.org/2000/svg">
+  <sodipodi:namedview
+     id="namedview346"
+     pagecolor="#ffffff"
+     bordercolor="#666666"
+     borderopacity="1.0"
+     inkscape:showpageshadow="2"
+     inkscape:pageopacity="0.0"
+     inkscape:pagecheckerboard="0"
+     inkscape:deskcolor="#d1d1d1"
+     inkscape:document-units="mm"
+     showgrid="false"
+     inkscape:zoom="5.6568543"
+     inkscape:cx="21.478368"
+     inkscape:cy="50.558134"
+     inkscape:window-width="1050"
+     inkscape:window-height="1619"
+     inkscape:window-x="0"
+     inkscape:window-y="25"
+     inkscape:window-maximized="1"
+     inkscape:current-layer="layer1" />
+  <defs
+     id="defs341">
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4450">
+      <g
+         id="g4456"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4454"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4452"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4466">
+      <g
+         id="g4472"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4470"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4468"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4474">
+      <g
+         id="g4480"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4478"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4476"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4482">
+      <g
+         id="g4488"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4486"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4484"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4490">
+      <g
+         id="g4496"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4494"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4492"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4498">
+      <g
+         id="g4504"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4502"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4500"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4514">
+      <g
+         id="g4520"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4518"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4516"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4530">
+      <g
+         id="g4536"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4534"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4532"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4546">
+      <g
+         id="g4552"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4550"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4548"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4554">
+      <g
+         id="g4560"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4558"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4556"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4562">
+      <g
+         id="g4568"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4566"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4564"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4570">
+      <g
+         id="g4576"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4574"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4572"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4578">
+      <g
+         id="g4584"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4582"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4580"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4586">
+      <g
+         id="g4592"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4590"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4588"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4594">
+      <g
+         id="g4600"
+         transform="matrix(0.03779528,0,0,-0.03779528,-1.0615633,3.4174527)">
+        <g
+           id="g4598"
+           transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+          <rect
+             x="30.000004"
+             y="-410"
+             width="360"
+             height="360"
+             id="rect4596"
+             style="stroke-width:100"
+             transform="scale(1,-1)" />
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4885">
+      <g
+         id="g4893"
+         transform="matrix(0.03929898,0,0,-0.02652666,-1.0787508,2.4265924)">
+        <g
+           id="g4891"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g4889"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect4887"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4895">
+      <g
+         id="g4903">
+        <g
+           id="g4901"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g4899"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect4897"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4925">
+      <g
+         id="g4933"
+         transform="matrix(0.03929898,0,0,-0.02652666,-1.0787508,2.4265924)">
+        <g
+           id="g4931"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g4929"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect4927"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4935">
+      <g
+         id="g4943">
+        <g
+           id="g4941"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g4939"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect4937"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4964">
+      <g
+         id="g4972"
+         transform="matrix(0.03929898,0,0,-0.02652666,-1.0787508,2.4265924)">
+        <g
+           id="g4970"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g4968"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect4966"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath4974">
+      <g
+         id="g4982">
+        <g
+           id="g4980"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g4978"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect4976"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <clipPath
+       clipPathUnits="userSpaceOnUse"
+       id="clipPath5120">
+      <g
+         id="g5128">
+        <g
+           id="g5126"
+           transform="matrix(0.96173703,0,0,1.4248035,0.43734998,-37.353372)">
+          <g
+             id="g5124"
+             transform="matrix(0.26458333,0,0,0.26458333,-27.475306,29.565937)">
+            <rect
+               x="30.000004"
+               y="-410"
+               width="360"
+               height="360"
+               id="rect5122"
+               style="stroke-width:100"
+               transform="scale(1,-1)" />
+          </g>
+        </g>
+      </g>
+    </clipPath>
+    <marker
+       style="overflow:visible"
+       refX="0"
+       refY="0"
+       orient="auto-start-reverse"
+       inkscape:stockid="Arrow2901564868668"
+       markerWidth="3.700125"
+       markerHeight="3.691"
+       viewBox="0 0 7.7 5.6"
+       inkscape:isstock="true"
+       inkscape:collect="always"
+       preserveAspectRatio="none"
+       id="Arrow2901564868668">
+      <path
+         transform="scale(0.7)"
+         d="M -2,-4 9,0 -2,4 c 2,-2.33 2,-5.66 0,-8 z"
+         style="fill:#ececec;fill-rule:evenodd;stroke:#ececec"
+         id="path6361" />
+    </marker>
+  </defs>
+  <g
+     inkscape:label="Layer 1"
+     inkscape:groupmode="layer"
+     id="layer1"
+     transform="translate(-71.260675,-78.317281)">
+    <rect
+       style="fill:#3d3d3d;fill-opacity:1;stroke:none;stroke-width:0;stroke-dasharray:none;stroke-opacity:1"
+       id="rect2922"
+       width="26.458328"
+       height="26.458334"
+       x="71.260681"
+       y="78.317284" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520831,0.18520831;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 84.511467,91.552827 c 0,0 1.035956,-1.077156 2.488832,-1.094515 1.452876,-0.01736 1.93503,0.569489 2.137195,0.813052 0.202166,0.243563 0.570948,1.09285 0.528038,1.640558 -0.04291,0.547709 -0.330294,1.455375 -0.758067,1.997564 -0.427774,0.542188 -1.077446,1.100601 -1.795551,1.325797 -0.718105,0.225196 -1.753121,0.161187 -2.26421,-0.31011 -0.511089,-0.471297 -0.993029,-1.327709 -0.939601,-2.105777 0.05343,-0.778069 0.06186,-1.020218 0.141182,-1.237757 0.07932,-0.217539 0.277603,-0.776094 0.462182,-1.028812 z"
+       id="path9078" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520797,0.18520797;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 84.483837,91.584443 c 0,0 1.841192,-1.914416 4.423368,-1.945267 2.582179,-0.03086 3.439106,1.012144 3.79841,1.445027 0.359306,0.432879 1.014739,1.942307 0.938477,2.915741 -0.07626,0.973436 -0.58703,2.586619 -1.347305,3.550245 -0.760278,0.963624 -1.914931,1.956084 -3.191209,2.356323 -1.276279,0.400238 -3.115798,0.286478 -4.024152,-0.551156 -0.908351,-0.83763 -1.764896,-2.35972 -1.669939,-3.742571 0.09497,-1.382851 0.109937,-1.81322 0.250922,-2.199848 0.14098,-0.38663 0.493379,-1.379343 0.821428,-1.828494 z"
+       id="path10328" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520831,0.18520831;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 84.546108,91.5411 c 0,0 -1.035956,1.077156 -2.488832,1.094515 -1.452876,0.01736 -1.93503,-0.569489 -2.137195,-0.813052 -0.202166,-0.243563 -0.570948,-1.09285 -0.528038,-1.640558 0.04291,-0.547709 0.330294,-1.455375 0.758067,-1.997564 0.427774,-0.542188 1.077446,-1.100601 1.795551,-1.325797 0.718105,-0.225196 1.753121,-0.161187 2.26421,0.31011 0.511089,0.471297 0.993029,1.327709 0.939601,2.105777 -0.05343,0.778069 -0.06186,1.020218 -0.141182,1.237757 -0.07932,0.217539 -0.277603,0.776094 -0.462182,1.028812 z"
+       id="path10330" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520797,0.18520797;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 84.573738,91.509484 c 0,0 -1.841192,1.914419 -4.423368,1.945269 -2.582179,0.0309 -3.439106,-1.012146 -3.79841,-1.445029 -0.359306,-0.432879 -1.014739,-1.942307 -0.938477,-2.915741 0.07626,-0.973436 0.58703,-2.586619 1.347305,-3.550245 0.760278,-0.963624 1.914931,-1.956084 3.191209,-2.356323 1.276279,-0.400237 3.115798,-0.286477 4.024152,0.551156 0.908351,0.83763 1.764896,2.35972 1.669939,3.742571 -0.09497,1.382851 -0.109937,1.81322 -0.250922,2.199848 -0.14098,0.38663 -0.493379,1.379343 -0.821428,1.828494 z"
+       id="path10332" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520831,0.18520831;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 94.584005,83.66225 c 0,0 -7.38303,4.522017 -10.154599,7.947337 -2.771568,3.425319 -5.600764,11.207533 -5.600764,11.207533"
+       id="path10334" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520831,0.18520831;stroke-dashoffset:0;stroke-opacity:1"
+       d="m 74.44805,99.353007 c 0,0 7.38303,-4.522015 10.154599,-7.947335 2.771568,-3.425319 5.600764,-11.207532 5.600764,-11.207532"
+       id="path10336" />
+    <path
+       style="fill:none;fill-opacity:1;stroke:#ababab;stroke-width:0.06173611;stroke-linecap:round;stroke-linejoin:miter;stroke-dasharray:0.18520831,0.18520831;stroke-dashoffset:0;stroke-opacity:1"
+       d="M 74.458074,103.8352 94.589302,79.288975"
+       id="path10338" />
+    <path
+       style="fill:#ececec;fill-opacity:1;stroke:#ececec;stroke-width:0.616;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1;marker-end:url(#Arrow2901564868668)"
+       d="m 81.837173,94.803608 4.389554,-5.365134"
+       id="path857"
+       sodipodi:nodetypes="cc" />
+    <circle
+       style="fill:#ececec;fill-opacity:1;stroke:none;stroke-width:0;stroke-dasharray:none;stroke-opacity:1"
+       id="circle3780"
+       cx="123.3643"
+       cy="17.329292"
+       r="0.98239028"
+       transform="rotate(39.28876)" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 76.419917,97.438267 0.599606,0.890159 0.559769,-1.024079 z"
+       id="path5218"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 78.280975,100.10051 -0.830674,-0.679627 1.071611,-0.462318 z"
+       id="path5214"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 80.175595,100.84119 -0.991276,-0.41145 0.892938,-0.751485 z"
+       id="path5216"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 80.212869,87.282831 0.830674,0.679627 -1.07161,0.462318 z"
+       id="path5212"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 76.866097,84.544627 0.830675,0.679627 -1.071611,0.462318 z"
+       id="path5210"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 89.562722,84.482341 0.830674,0.679627 0.240937,-1.141944 z"
+       id="path5220"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 91.557688,85.002596 0.599607,0.89016 0.559769,-1.02408 z"
+       id="path5222"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 89.523835,83.332525 -0.991276,-0.411441 0.892938,-0.751495 z"
+       id="path5224"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 88.716453,94.240149 0.830674,0.679627 -1.071611,0.462318 z"
+       id="path5206"
+       sodipodi:nodetypes="cccc" />
+    <path
+       style="fill:#ababab;fill-opacity:1;stroke:none;stroke-width:0.0445293;stroke-linecap:round;stroke-dasharray:none;stroke-opacity:1"
+       d="m 92.03266,96.953347 0.830674,0.679627 -1.07161,0.462318 z"
+       id="path5208"
+       sodipodi:nodetypes="cccc" />
+  </g>
+</svg>
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 6a00a12..860fe90 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -65,12 +65,26 @@ Equivalent Sources
 Forward modelling
 -----------------
 
+Gravity fields:
+
 .. autosummary::
     :toctree: generated/
 
     point_gravity
     prism_gravity
     tesseroid_gravity
+
+Magnetic fields:
+
+.. autosummary::
+    :toctree: generated/
+
+    dipole_magnetic
+    dipole_magnetic_component
+
+Layers and meshes:
+
+.. autosummary::
     prism_layer
     tesseroid_layer
     DatasetAccessorPrismLayer
diff --git a/doc/user_guide/forward_modelling/dipole.rst b/doc/user_guide/forward_modelling/dipole.rst
new file mode 100644
index 0000000..e4f0dee
--- /dev/null
+++ b/doc/user_guide/forward_modelling/dipole.rst
@@ -0,0 +1,123 @@
+.. _dipole:
+
+Magnetic dipoles
+================
+
+We can compute the magnetic field generated by dipole sources through
+:func:`harmonica.dipole_magnetic` and
+:func:`harmonica.dipole_magnetic_component` functions.
+The former computes the three components of the magnetic field generated by
+a set of dipoles on every observation point, while the latter computes only
+a single component.
+
+Each dipole can be defined as a tuple containing its coordinates in the
+following order: *easting*, *northing*, *upward* (in Cartesian coordinates and
+in meters).
+
+.. jupyter-execute::
+   :hide-code:
+
+   import harmonica as hm
+
+Let's define a single dipole and compute the magnetic field it generates on
+a computation point above it.
+
+Define the dipole and its magnetic moment vector:
+
+.. jupyter-execute::
+
+   dipole = (20, 40, -50)
+   magnetic_moment = [100, 100, 100]
+
+Define the observation point above it:
+
+.. jupyter-execute::
+
+   coordinates = (20, 40, 10)
+
+And compute the three components of the magnetic field the dipole generates:
+
+.. jupyter-execute::
+
+   b_e, b_n, b_u = hm.dipole_magnetic(coordinates, dipole, magnetic_moment)
+   print(b_e, b_n, b_u)
+
+We can use :func:`harmonica.dipole_magnetic_component` to compute just a single
+component:
+
+.. jupyter-execute::
+
+   b_u = hm.dipole_magnetic_component(
+      coordinates, dipole, magnetic_moment, component="upward"
+   )
+   print(b_u)
+
+
+Multiple dipoles
+----------------
+
+These two functions also allow us to compute the magnetic fields of multiple
+dipoles on multiple observation points.
+
+Let's define a regular grid of observation points:
+
+.. jupyter-execute::
+
+   import verde as vd
+
+   region = (-100, 100, -100, 100)
+   spacing = 1
+   height = 0
+   coordinates = vd.grid_coordinates(
+      region=region, spacing=spacing, extra_coords=height
+   )
+
+And a set of dipoles with their magnetic moments:
+
+.. jupyter-execute::
+
+   import numpy as np
+
+   easting = [25, 35, -30, -50]
+   northing = [3, -38, 22, -30]
+   upward = [-200, -100, -300, -150]
+   dipoles = (easting, northing, upward)
+
+   magnetic_moments = np.array([
+      [1e3, 1e3, 1e3],
+      [2e3, 2e3, 2e3],
+      [500, 500, 500],
+      [2e3, 2e3, 2e3],
+   ])
+
+Now, let's compute the magnetic field components that the dipoles generate on
+every observation point:
+
+.. jupyter-execute::
+
+   b_e, b_n, b_u = hm.dipole_magnetic(coordinates, dipoles, magnetic_moments)
+
+.. jupyter-execute::
+
+   import matplotlib.pyplot as plt
+
+   fig, axes = plt.subplots(nrows=1, ncols=3, sharey=True, figsize=(12, 8))
+
+   fields = {"b_e": b_e, "b_n": b_n, "b_u": b_u}
+   for field, ax in zip(fields, axes):
+      tmp = ax.pcolormesh(coordinates[0], coordinates[1], fields[field])
+      ax.set_aspect("equal")
+      ax.set_title(field)
+      ax.ticklabel_format(style="sci", scilimits=(0, 0), axis="both")
+      plt.colorbar(tmp, ax=ax, orientation="horizontal", label="nT", pad=0.008)
+   plt.show()
+
+----
+
+.. grid:: 2
+
+    .. grid-item-card:: :jupyter-download-script:`Download Python script <dipole>`
+        :text-align: center
+
+    .. grid-item-card:: :jupyter-download-nb:`Download Jupyter notebook <dipole>`
+        :text-align: center
diff --git a/doc/user_guide/forward_modelling/index.rst b/doc/user_guide/forward_modelling/index.rst
index cbd023e..554136f 100644
--- a/doc/user_guide/forward_modelling/index.rst
+++ b/doc/user_guide/forward_modelling/index.rst
@@ -8,16 +8,15 @@ Forward modelling
    :hidden:
 
    point
+   dipole
    prism
    tesseroid
 
 Harmonica offers some functions and classes to calculate the gravity and
 magnetic fields generated by the following 3D bodies.
 
-Gravity
--------
-
-.. grid:: 1 1 3 3
+.. grid:: 1 1 2 2
+   :gutter: 2 3 5 5
 
    .. grid-item-card:: Point sources
        :img-top: ../../_static/figures/simple_points.svg
@@ -27,6 +26,14 @@ Gravity
 
        Cartesian and spherical coordinates
 
+   .. grid-item-card:: Magnetic dipoles
+       :img-top: ../../_static/figures/dipoles.svg
+       :text-align: center
+       :link: dipole
+       :link-type: ref
+
+       Cartesian coordinates
+
    .. grid-item-card:: Rectangular prisms
        :img-top: ../../_static/figures/simple_prism.svg
        :text-align: center
diff --git a/harmonica/__init__.py b/harmonica/__init__.py
index 8bd122a..33dc1b2 100644
--- a/harmonica/__init__.py
+++ b/harmonica/__init__.py
@@ -9,6 +9,7 @@
 from ._equivalent_sources.cartesian import EquivalentSources
 from ._equivalent_sources.gradient_boosted import EquivalentSourcesGB
 from ._equivalent_sources.spherical import EquivalentSourcesSph
+from ._forward.dipole import dipole_magnetic, dipole_magnetic_component
 from ._forward.point import point_gravity
 from ._forward.prism import prism_gravity
 from ._forward.prism_layer import DatasetAccessorPrismLayer, prism_layer
diff --git a/harmonica/_forward/dipole.py b/harmonica/_forward/dipole.py
new file mode 100644
index 0000000..b3567e5
--- /dev/null
+++ b/harmonica/_forward/dipole.py
@@ -0,0 +1,309 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Forward modelling for magnetic fields of dipoles
+"""
+import numpy as np
+from choclo.dipole import magnetic_e, magnetic_field, magnetic_n, magnetic_u
+from numba import jit, prange
+
+
+def dipole_magnetic(
+    coordinates,
+    dipoles,
+    magnetic_moments,
+    parallel=True,
+    dtype="float64",
+    disable_checks=False,
+):
+    """
+    Magnetic field of dipoles in Cartesian coordinates
+
+    Parameters
+    ----------
+    coordinates : list of arrays
+        List of arrays containing the ``easting``, ``northing`` and ``upward``
+        coordinates of the computation points defined on a Cartesian coordinate
+        system. All coordinates should be in meters.
+    dipoles : list of arrays
+        List of arrays containing the ``easting``, ``northing`` and ``upward``
+        locations of the dipoles defined on a Cartesian coordinate system. All
+        coordinates should be in meters.
+    magnetic_moments : 2d-array
+        Two dimensional array containing the magnetic moments of each dipole in
+        :math:`Am^2`. Each row contains the components of the magnetic moment
+        for each dipole in the following order: ``mag_moment_easting``,
+        ``mag_moment_northing``, ``mag_moment_upward``.
+    parallel : bool (optional)
+        If True the computations will run in parallel using Numba built-in
+        parallelization. If False, the forward model will run on a single core.
+        Might be useful to disable parallelization if the forward model is run
+        by an already parallelized workflow. Default to True.
+    dtype : data-type (optional)
+        Data type assigned to the resulting magnetic field. Default to
+        ``np.float64``.
+    disable_checks : bool (optional)
+        Flag that controls whether to perform a sanity check on the model.
+        Should be set to ``True`` only when it is certain that the input model
+        is valid and it does not need to be checked.
+        Default to ``False``.
+
+    Returns
+    -------
+    magnetic_field : tuple of arrays
+        Tuple containing each component of the magnetic field generated by the
+        dipoles as arrays. The three components are returned in the following
+        order: ``b_e``, ``b_n``, ``b_u``. Magnetic field components are in
+        :math:`nT`.
+    """
+    # Figure out the shape and size of the output arrays
+    cast = np.broadcast(*coordinates[:3])
+    b_e, b_n, b_u = tuple(np.zeros(cast.size, dtype=dtype) for _ in range(3))
+    # Convert coordinates, dipoles and magnetic moments to arrays
+    coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3])
+    dipoles = tuple(np.atleast_1d(i).ravel() for i in dipoles[:3])
+    magnetic_moments = np.atleast_2d(magnetic_moments)
+    # Sanity checks
+    if not disable_checks:
+        _check_dipoles_and_magnetic_moments(dipoles, magnetic_moments)
+    # Compute the magnetic fields of the dipoles
+    if parallel:
+        _jit_dipole_magnetic_field_cartesian_parallel(
+            coordinates, dipoles, magnetic_moments, b_e, b_n, b_u
+        )
+    else:
+        _jit_dipole_magnetic_field_cartesian_serial(
+            coordinates, dipoles, magnetic_moments, b_e, b_n, b_u
+        )
+    # Convert to nT
+    b_e *= 1e9
+    b_n *= 1e9
+    b_u *= 1e9
+    return b_e.reshape(cast.shape), b_n.reshape(cast.shape), b_u.reshape(cast.shape)
+
+
+def dipole_magnetic_component(
+    coordinates,
+    dipoles,
+    magnetic_moments,
+    component,
+    parallel=True,
+    dtype="float64",
+    disable_checks=False,
+):
+    """
+    Single magnetic field component of dipoles in Cartesian coordinates
+
+    .. important::
+
+        Use this function only if you need to compute a single component of the
+        magnetic field. Use :func:`harmonica.dipole_magnetic` to compute the
+        three components more efficiently.
+
+    Parameters
+    ----------
+    coordinates : list of arrays
+        List of arrays containing the ``easting``, ``northing`` and ``upward``
+        coordinates of the computation points defined on a Cartesian coordinate
+        system. All coordinates should be in meters.
+    dipoles : list of arrays
+        List of arrays containing the ``easting``, ``northing`` and ``upward``
+        locations of the dipoles defined on a Cartesian coordinate system. All
+        coordinates should be in meters.
+    magnetic_moments : 2d-array
+        Two dimensional array containing the magnetic moments of each dipole in
+        :math:`Am^2`. Each row contains the components of the magnetic moment
+        for each dipole in the following order: ``mag_moment_easting``,
+        ``mag_moment_northing``, ``mag_moment_upward``.
+    component : str
+        Computed that will be computed. Available options are: ``"easting"``,
+        ``"northing"`` or ``"upward"``.
+    parallel : bool (optional)
+        If True the computations will run in parallel using Numba built-in
+        parallelization. If False, the forward model will run on a single core.
+        Might be useful to disable parallelization if the forward model is run
+        by an already parallelized workflow. Default to True.
+    dtype : data-type (optional)
+        Data type assigned to the resulting magnetic field. Default to
+        ``np.float64``.
+    disable_checks : bool (optional)
+        Flag that controls whether to perform a sanity check on the model.
+        Should be set to ``True`` only when it is certain that the input model
+        is valid and it does not need to be checked.
+        Default to ``False``.
+
+    Returns
+    -------
+    magnetic_field :  array
+        Array containing the magnetic field component in :math:`nT`.
+    """
+    # Figure out the shape and size of the output arrays
+    cast = np.broadcast(*coordinates[:3])
+    result = np.zeros(cast.size, dtype=dtype)
+    # Convert coordinates, dipoles and magnetic moments to arrays
+    coordinates = tuple(np.atleast_1d(i).ravel() for i in coordinates[:3])
+    dipoles = tuple(np.atleast_1d(i).ravel() for i in dipoles[:3])
+    magnetic_moments = np.atleast_2d(magnetic_moments)
+    # Choose forward modelling function based on the chosen component
+    forward_function = _get_magnetic_forward_function(component)
+    # Sanity checks
+    if not disable_checks:
+        _check_dipoles_and_magnetic_moments(dipoles, magnetic_moments)
+    # Compute the magnetic fields of the dipoles
+    if parallel:
+        _jit_dipole_magnetic_component_cartesian_parallel(
+            coordinates, dipoles, magnetic_moments, result, forward_function
+        )
+    else:
+        _jit_dipole_magnetic_component_cartesian_serial(
+            coordinates, dipoles, magnetic_moments, result, forward_function
+        )
+    # Convert to nT
+    result *= 1e9
+    return result.reshape(cast.shape)
+
+
+def _check_dipoles_and_magnetic_moments(dipoles, magnetic_moments):
+    """
+    Check if dipoles and magnetic moments have valid shape and size
+    """
+    if magnetic_moments.shape[1] != 3:
+        raise ValueError(
+            f"Invalid magnetic moments with '{magnetic_moments.shape[1]}' elements."
+            " Magnetic moments vectors should have 3 components."
+        )
+    if magnetic_moments.shape[0] != dipoles[0].size:
+        raise ValueError(
+            f"Number of elements in magnetic_moments ({magnetic_moments.shape[0]})"
+            f" mismatch the number of dipoles ({dipoles[0].size})."
+        )
+
+
+def _jit_dipole_magnetic_field_cartesian(
+    coordinates, dipoles, magnetic_moments, b_e, b_n, b_u
+):
+    """
+    Compute the magnetic field components of dipoles in Cartesian coordinates
+
+    Parameters
+    ----------
+    coordinates : tuple
+        Tuple containing the ``easting``, ``northing``, ``upward`` coordinates
+        of the computation points in Cartesian coordinate system.
+    dipoles : tuple
+        Tuple containing the ``easting``, ``northing``, ``upward`` coordinates
+        of the dipoles in Cartesian coordinate system.
+    magnetic_moments : 2d-array
+        Array with the components of the magnetic moment for each dipole.
+    b_e : 1d-array
+        Array where the resulting values of the easting component of the
+        magnetic field will be stored.
+    b_n : 1d-array
+        Array where the resulting values of the northing component of the
+        magnetic field will be stored.
+    b_u : 1d-array
+        Array where the resulting values of the upward component of the
+        magnetic field will be stored.
+    """
+    easting, northing, upward = coordinates
+    easting_p, northing_p, upward_p = dipoles
+    for l in prange(easting.size):
+        for m in range(easting_p.size):
+            easting_comp, northing_comp, upward_comp = magnetic_field(
+                easting[l],
+                northing[l],
+                upward[l],
+                easting_p[m],
+                northing_p[m],
+                upward_p[m],
+                magnetic_moments[m, 0],
+                magnetic_moments[m, 1],
+                magnetic_moments[m, 2],
+            )
+            b_e[l] += easting_comp
+            b_n[l] += northing_comp
+            b_u[l] += upward_comp
+
+
+def _jit_dipole_magnetic_component_cartesian(
+    coordinates, dipoles, magnetic_moments, result, forward_func
+):
+    """
+    Compute a single magnetic component of dipoles in Cartesian coordinates
+
+    Parameters
+    ----------
+    coordinates : tuple
+        Tuple containing the ``easting``, ``northing``, ``upward`` coordinates
+        of the computation points in Cartesian coordinate system.
+    dipoles : tuple
+        Tuple containing the ``easting``, ``northing``, ``upward`` coordinates
+        of the dipoles in Cartesian coordinate system.
+    magnetic_moments : 2d-array
+        Array with the components of the magnetic moment for each dipole.
+    result : 1d-array
+        Array where the resulting values of the selected component of the
+        magnetic field will be stored.
+    forward_function : callable
+        Forward function to be used to compute the desired component of the
+        magnetic field. Choose one of :func:`choclo.dipole.magnetic_easting`,
+        :func:`choclo.dipole.magnetic_northing` or
+        :func:`choclo.dipole.magnetic_upward`.
+    """
+    easting, northing, upward = coordinates
+    easting_p, northing_p, upward_p = dipoles
+    for l in prange(easting.size):
+        for m in range(easting_p.size):
+            result[l] += forward_func(
+                easting[l],
+                northing[l],
+                upward[l],
+                easting_p[m],
+                northing_p[m],
+                upward_p[m],
+                magnetic_moments[m, 0],
+                magnetic_moments[m, 1],
+                magnetic_moments[m, 2],
+            )
+
+
+def _get_magnetic_forward_function(component):
+    """
+    Returns the Choclo magnetic forward modelling function for the desired
+    component
+
+    Parameters
+    ----------
+    component : str
+        Magnetic field component.
+
+    Returns
+    -------
+    forward_function : callable
+        Forward modelling function for the desired component.
+    """
+    if component not in ("easting", "northing", "upward"):
+        raise ValueError(
+            f"Invalid component '{component}'. "
+            "It must be either 'easting', 'northing' or 'upward'."
+        )
+    functions = {"easting": magnetic_e, "northing": magnetic_n, "upward": magnetic_u}
+    return functions[component]
+
+
+_jit_dipole_magnetic_field_cartesian_serial = jit(nopython=True)(
+    _jit_dipole_magnetic_field_cartesian
+)
+_jit_dipole_magnetic_field_cartesian_parallel = jit(nopython=True, parallel=True)(
+    _jit_dipole_magnetic_field_cartesian
+)
+_jit_dipole_magnetic_component_cartesian_serial = jit(nopython=True)(
+    _jit_dipole_magnetic_component_cartesian
+)
+_jit_dipole_magnetic_component_cartesian_parallel = jit(nopython=True, parallel=True)(
+    _jit_dipole_magnetic_component_cartesian
+)


Test patch:
diff --git a/harmonica/tests/test_dipole.py b/harmonica/tests/test_dipole.py
new file mode 100644
index 0000000..8549518
--- /dev/null
+++ b/harmonica/tests/test_dipole.py
@@ -0,0 +1,284 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Test magnetic forward functions for dipoles
+"""
+import choclo
+import numpy as np
+import numpy.testing as npt
+import pytest
+import verde as vd
+from choclo.dipole import magnetic_field
+
+from .. import dipole_magnetic, dipole_magnetic_component
+from .utils import run_only_with_numba
+
+
+def test_invalid_component():
+    "Check if passing an invalid component raises an error"
+    coordinates = [0, 0, 0]
+    dipoles = [[-100, 100], [-100, 100], [-200, -100]]
+    magnetic_moments = [[1, 1, 2], [-1, -2, 3]]
+    with pytest.raises(ValueError, match="Invalid component"):
+        dipole_magnetic_component(
+            coordinates, dipoles, magnetic_moments, component="Not a valid field"
+        )
+
+
+class TestInvalidMagneticMoments:
+    @pytest.mark.parametrize("component", ("all", "easting", "northing", "upward"))
+    def test_magnetic_moments_and_dipoles(self, component):
+        """
+        Test error when dipoles and magnetic moments mismatch
+        """
+        coordinates = [0, 0, 0]
+        dipoles = [[1, 1], [-1, 0], [-10, -2]]
+        # Define invalid magnetic moments (more moments than dipoles)
+        magnetic_moments = [
+            [1, 2, 3],
+            [4, 5, 6],
+            [7, 8, 9],
+        ]
+        msg = "Number of elements in magnetic_moments"
+        if component == "all":
+            with pytest.raises(ValueError, match=msg):
+                dipole_magnetic(coordinates, dipoles, magnetic_moments)
+        else:
+            with pytest.raises(ValueError, match=msg):
+                dipole_magnetic_component(
+                    coordinates, dipoles, magnetic_moments, component
+                )
+
+    @pytest.mark.parametrize("component", ("all", "easting", "northing", "upward"))
+    def test_magnetic_moments_components(self, component):
+        """
+        Test if error is raised when magnetic moments have != 3 components
+        """
+        coordinates = [0, 0, 0]
+        dipoles = [[1, 1], [-1, 0], [-10, -2]]
+        # Define invalid magnetic moments (more than three components)
+        magnetic_moments = [
+            [1, 2, 3, 5],
+            [4, 5, 6, 7],
+        ]
+        msg = "Invalid magnetic moments with '4' elements."
+        if component == "all":
+            with pytest.raises(ValueError, match=msg):
+                dipole_magnetic(coordinates, dipoles, magnetic_moments)
+        else:
+            with pytest.raises(ValueError, match=msg):
+                dipole_magnetic_component(
+                    coordinates, dipoles, magnetic_moments, component
+                )
+
+    @pytest.mark.parametrize("component", ("all", "easting", "northing", "upward"))
+    def test_disable_checks(self, component):
+        """
+        Test disable checks on invalid magnetic moments
+        """
+        coordinates = [0, 0, 0]
+        dipoles = [[1, 1], [-1, 0], [-10, -2]]
+        # Define invalid magnetic moments (more than three components)
+        magnetic_moments = [
+            [1, 2, 3, 5],
+            [4, 5, 6, 7],
+        ]
+        if component == "all":
+            dipole_magnetic(coordinates, dipoles, magnetic_moments, disable_checks=True)
+        else:
+            dipole_magnetic_component(
+                coordinates,
+                dipoles,
+                magnetic_moments,
+                component,
+                disable_checks=True,
+            )
+
+
+class TestSerialVsParallel:
+    """
+    Test serial vs parallel
+    """
+
+    @pytest.mark.parametrize("component", ("all", "easting", "northing", "upward"))
+    def test_dipoles_parallel_vs_serial_no_numba(self, component):
+        """
+        Check results of parallelized and serials runs
+        Run a small problem with Numba disable to count for test coverage.
+        """
+        coordinates = ([0, 10], [0, 10], [0, 10])
+        dipoles = [[-100, 0], [0, 100], [-20, -50]]
+        magnetic_moments = [
+            [1.0, 1.0, 1.0],
+            [1.0, -1.0, 5.0],
+        ]
+        if component == "all":
+            parallel = dipole_magnetic(
+                coordinates, dipoles, magnetic_moments, parallel=True
+            )
+            serial = dipole_magnetic(
+                coordinates, dipoles, magnetic_moments, parallel=False
+            )
+        else:
+            parallel = dipole_magnetic_component(
+                coordinates, dipoles, magnetic_moments, component, parallel=True
+            )
+            serial = dipole_magnetic_component(
+                coordinates, dipoles, magnetic_moments, component, parallel=False
+            )
+        npt.assert_allclose(parallel, serial)
+
+    @run_only_with_numba
+    @pytest.mark.parametrize("component", ("all", "easting", "northing", "upward"))
+    def test_dipoles_parallel_vs_serial(self, component):
+        """
+        Check results of parallelized and serials runs
+        Run a large problem only with Numba enabled.
+        """
+        coordinates = vd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, extra_coords=10
+        )
+        dipoles = [[-100, 0], [0, 100], [-20, -50]]
+        magnetic_moments = [
+            [1.0, 1.0, 1.0],
+            [1.0, -1.0, 5.0],
+        ]
+        if component == "all":
+            parallel = dipole_magnetic(
+                coordinates, dipoles, magnetic_moments, parallel=True
+            )
+            serial = dipole_magnetic(
+                coordinates, dipoles, magnetic_moments, parallel=False
+            )
+        else:
+            parallel = dipole_magnetic_component(
+                coordinates, dipoles, magnetic_moments, component, parallel=True
+            )
+            serial = dipole_magnetic_component(
+                coordinates, dipoles, magnetic_moments, component, parallel=False
+            )
+        npt.assert_allclose(parallel, serial)
+
+
+class TestAgainstChoclo:
+    """
+    Test forward modelling functions against dumb Choclo runs
+    """
+
+    @pytest.fixture()
+    def sample_dipoles(self):
+        """
+        Return four sample dipoles
+        """
+        dipoles = (
+            np.array([-10.0, 0.0, 10.0, 5.0]),
+            np.array([2.0, 0.0, -1.0, 1.5]),
+            np.array([-10.0, -15.0, -5.0, -2.0]),
+        )
+        return dipoles
+
+    @pytest.fixture()
+    def sample_magnetic_moments(self):
+        """
+        Return sample magnetic moment vectors for the dipoles
+        """
+        magnetic_moments = np.array(
+            [
+                [1.0, 1.0, 1],
+                [-1.0, -0.5, 2.0],
+                [3, -1.0, -3.0],
+                [-1.5, 2.0, -2.0],
+            ],
+            dtype=np.float64,
+        )
+        return magnetic_moments
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return four sample observation points
+        """
+        easting = np.array([-5, 10, 0, 15], dtype=float)
+        northing = np.array([14, -4, 11, 0], dtype=float)
+        upward = np.array([9, 6, 6, 12], dtype=float)
+        return (easting, northing, upward)
+
+    def test_against_choclo(
+        self, sample_coordinates, sample_dipoles, sample_magnetic_moments
+    ):
+        """
+        Test dipole_magnetic against raw Choclo runs
+        """
+        easting, northing, upward = sample_coordinates
+        easting_p, northing_p, upward_p = sample_dipoles
+        # Compute expected results with dumb Choclo runs
+        n_coords = easting.size
+        n_dipoles = sample_dipoles[0].size
+        expected_magnetic_e = np.zeros(n_coords, dtype=np.float64)
+        expected_magnetic_n = np.zeros(n_coords, dtype=np.float64)
+        expected_magnetic_u = np.zeros(n_coords, dtype=np.float64)
+        for i in range(n_coords):
+            for j in range(n_dipoles):
+                b_e, b_n, b_u = magnetic_field(
+                    easting[i],
+                    northing[i],
+                    upward[i],
+                    easting_p[j],
+                    northing_p[j],
+                    upward_p[j],
+                    *sample_magnetic_moments[j, :],
+                )
+                expected_magnetic_e[i] += b_e
+                expected_magnetic_n[i] += b_n
+                expected_magnetic_u[i] += b_u
+        # Convert to nT
+        expected_magnetic_e *= 1e9
+        expected_magnetic_n *= 1e9
+        expected_magnetic_u *= 1e9
+        # Compare with harmonica results
+        b_e, b_n, b_u = dipole_magnetic(
+            sample_coordinates, sample_dipoles, sample_magnetic_moments
+        )
+        npt.assert_allclose(b_e, expected_magnetic_e)
+        npt.assert_allclose(b_n, expected_magnetic_n)
+        npt.assert_allclose(b_u, expected_magnetic_u)
+
+    @pytest.mark.parametrize("component", ("easting", "northing", "upward"))
+    def test_component_against_choclo(
+        self, sample_coordinates, sample_dipoles, sample_magnetic_moments, component
+    ):
+        """
+        Test dipole_magnetic_component against raw Choclo runs
+        """
+        easting, northing, upward = sample_coordinates
+        easting_p, northing_p, upward_p = sample_dipoles
+        # Compute expected results with dumb Choclo runs
+        n_coords = easting.size
+        n_dipoles = sample_dipoles[0].size
+        expected_result = np.zeros(n_coords, dtype=np.float64)
+        forward_func = getattr(choclo.dipole, f"magnetic_{component[0]}")
+        for i in range(n_coords):
+            for j in range(n_dipoles):
+                expected_result[i] += forward_func(
+                    easting[i],
+                    northing[i],
+                    upward[i],
+                    easting_p[j],
+                    northing_p[j],
+                    upward_p[j],
+                    *sample_magnetic_moments[j, :],
+                )
+        # Convert to nT
+        expected_result *= 1e9
+        # Compare with harmonica results
+        result = dipole_magnetic_component(
+            sample_coordinates,
+            sample_dipoles,
+            sample_magnetic_moments,
+            component,
+        )
+        npt.assert_allclose(result, expected_result)