diff --git a/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourcePlanar.java b/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourcePlanar.java
new file mode 100644
index 0000000000000000000000000000000000000000..6d82e218598e231f73fd60f601007c47c5599856
--- /dev/null
+++ b/src/main/java/gov/usgs/earthquake/nshmp/model/GridSourcePlanar.java
@@ -0,0 +1,288 @@
+package gov.usgs.earthquake.nshmp.model;
+
+import static gov.usgs.earthquake.nshmp.Maths.PI_BY_2;
+import static gov.usgs.earthquake.nshmp.Maths.TO_RADIANS;
+import static gov.usgs.earthquake.nshmp.Maths.TWO_PI;
+import static gov.usgs.earthquake.nshmp.Maths.hypot;
+import static java.lang.Math.min;
+import static java.lang.Math.sin;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.tan;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;
+import java.util.Optional;
+
+import com.google.common.math.DoubleMath;
+
+import gov.usgs.earthquake.nshmp.data.XySequence;
+import gov.usgs.earthquake.nshmp.fault.FocalMech;
+import gov.usgs.earthquake.nshmp.fault.surface.RuptureScaling;
+import gov.usgs.earthquake.nshmp.fault.surface.RuptureSurface;
+import gov.usgs.earthquake.nshmp.geo.Location;
+import gov.usgs.earthquake.nshmp.geo.LocationVector;
+import gov.usgs.earthquake.nshmp.geo.Locations;
+import gov.usgs.earthquake.nshmp.mfd.Mfd;
+import gov.usgs.earthquake.nshmp.model.GridLoader.RuptureData;
+
+/**
+ * Planar gridded seismicity source. For initial use with subduction gridded
+ * seismicity.
+ *
+ * The internal MFD will go up to Mmax for the entire parent rupture set, but
+ * the highest magnitudes will often have zero rate.
+ *
+ * @author U.S. Geological Survey
+ */
+class GridSourcePlanar extends PointSource {
+
+ final RuptureSet<? extends Source> parent;
+ final SourceType type;
+ final Location loc;
+ final Mfd mfd;
+ final List<Rupture> ruptures;
+
+ /**
+ * Constructs a new planar grid earthquake source. This is a simple model that
+ * does not simulate finiteness (e.g. rupture rRup values will differ from rJB
+ * only by virtue of the depth of the source).
+ *
+ * @param parent {@code RuptureSet}
+ * @param type of source, as supplied from a parent {@code RuptureSet}
+ * @param loc <code>Location</code> of the point source
+ * @param mfd magnitude frequency distribution of the source
+ * @param mechWtMap <code>Map</code> of focal mechanism weights
+ * @param rupScaling rupture scaling model that may, or may not, impose an rJB
+ * distance correction
+ * @param depthModel specifies magnitude cutoffs and associated weights for
+ * different depth-to-top-of-ruptures
+ * @param surface optional rupture surface instance for use by subclasses
+ */
+ GridSourcePlanar(
+ RuptureSet<PointSource> parent,
+ SourceType type,
+ Location loc,
+ Mfd mfd,
+ Map<FocalMech, Double> mechWtMap,
+ RuptureScaling rupScaling,
+ DepthModel depthModel,
+ List<RuptureData> ruptureData) {
+
+ super(parent, type, loc, mfd, mechWtMap, rupScaling, depthModel, Optional.empty());
+
+ this.parent = parent;
+ this.type = type;
+ this.loc = loc;
+ this.mfd = mfd;
+
+ this.ruptures = initRuptures(loc, mfd, ruptureData);
+ }
+
+ private List<Rupture> initRuptures(
+ Location loc,
+ Mfd mfd,
+ List<RuptureData> ruptureData) {
+
+ List<Rupture> ruptures = new ArrayList<>(ruptureData.size());
+ XySequence mr = mfd.data();
+ for (int i = 0; i < ruptureData.size(); i++) {
+ PlanarSurface surface = new PlanarSurface(loc, ruptureData.get(i));
+ Rupture rupture = new RegularRupture(mr.x(i), mr.y(i), 0, surface);
+ ruptures.add(rupture);
+ }
+ return ruptures;
+ }
+
+ @Override
+ public String name() {
+ return "GridSourcePlanar: " + formatLocation(loc);
+ }
+
+ static final String FORMAT = "%.3f, %.3f";
+
+ static String formatLocation(Location loc) {
+ return String.format(FORMAT, loc.longitude, loc.latitude);
+ }
+
+ @Override
+ public int size() {
+ return ruptures.size();
+ }
+
+ @Override
+ public Rupture get(int index) {
+ return ruptures.get(index);
+ }
+
+ /**
+ * Overridden to return the ID of the parent grid rupture set. All point
+ * sources instantiated for a grid source will have the same ID. Note that
+ * {@code PointSource}s may be retrieved by index using
+ * {@link GridRuptureSet#source(int)}.
+ */
+ @Override
+ public int id() {
+ return parent.id();
+ }
+
+ @Override
+ public SourceType type() {
+ return type;
+ }
+
+ /**
+ * The location of the point source, irrespective of any distance corrections
+ * that might be applied to attendant ruptures and ignoring the supplied site
+ * {@code Location}.
+ */
+ @Override
+ public Location location(Location site) {
+ return loc;
+ }
+
+ @Override
+ public List<Mfd> mfds() {
+ return List.of(mfd);
+ }
+
+ static class PlanarSurface implements RuptureSurface {
+
+ final Location location;
+ final double strike;
+ final double dip;
+ final double dipRad;
+ final double zTor;
+ final double zBor;
+ final double length;
+ final double width;
+ final double widthH;
+
+ // Four corners of rupture:
+ // top trace: p1 --> p2
+ // bot trace: p4 <-- p3
+ Location p1, p2, p3, p4;
+
+ PlanarSurface(Location location, RuptureData rupData) {
+
+ this.location = location;
+
+ this.strike = rupData.strike;
+ this.dip = rupData.dip;
+ this.dipRad = dip * TO_RADIANS;
+ this.zTor = rupData.zTor;
+ this.zBor = rupData.zBor;
+ this.length = rupData.length;
+
+ double strikeRad = strike * TO_RADIANS;
+
+ // forward and back strike vectors
+ LocationVector svF = LocationVector.create(strikeRad, length / 2.0, 0.0);
+ LocationVector svB = LocationVector.reverseOf(svF);
+
+ // down-dip and up-dip azimuths
+ double ddAz = (strikeRad + PI_BY_2) % TWO_PI;
+ double udAz = (strikeRad - PI_BY_2) % TWO_PI;
+
+ // down-dip
+ double ddΔV = zBor - location.depth;
+ double ddΔH = ddΔV / tan(dipRad);
+ double ddW = ddΔV / sin(dipRad);
+ LocationVector ddV = LocationVector.create(ddAz, ddΔH, -ddΔV);
+ Location ddP = Locations.location(location, ddV);
+
+ // up-dip
+ double udΔV = location.depth - zTor;
+ double udΔH = udΔV / tan(dipRad);
+ double udW = udΔV / sin(dipRad);
+ LocationVector udV = LocationVector.create(udAz, udΔH, udΔV);
+ Location udP = Locations.location(location, udV);
+
+ this.p1 = Locations.location(udP, svB);
+ this.p2 = Locations.location(udP, svF);
+ this.p3 = Locations.location(ddP, svF);
+ this.p4 = Locations.location(ddP, svB);
+
+ this.width = ddW + udW;
+ this.widthH = ddΔH + udΔH;
+ }
+
+ @Override
+ public Distance distanceTo(Location loc) {
+ // NOTE no NSHMP style distance corrections here
+
+ double rX = Locations.distanceToLineFast(p1, p2, loc);
+ double rSeg = Locations.distanceToSegmentFast(p1, p2, loc);
+
+ // simple footwall case
+ boolean isVertical = (dip() > 89.0);
+ if (rX <= 0.0 || isVertical) {
+ return Distance.create(rSeg, hypot(rSeg, zTor), rX);
+ }
+
+ // otherwise, we're on the hanging wall...
+
+ // compute rRup as though we're between endpoints
+ double rCutTop = tan(dipRad) * zTor;
+ double rCutBot = tan(dipRad) * zBor + widthH;
+ double rRup = (rX > rCutBot) ? hypot(rX - widthH, zBor) : (rX < rCutTop) ? hypot(rX,
+ zTor) : hypot(rCutTop, zTor) + (rX - rCutTop) * sin(dipRad);
+
+ // test if we're normal to trace or past its endpoints
+ boolean offEnd = DoubleMath.fuzzyCompare(rSeg, rX, 0.00001) > 0;
+
+ if (offEnd) {
+ // distance from surface projection of ends/caps of fault
+ double rJB = min(Locations.distanceToSegmentFast(p1, p4, loc),
+ Locations.distanceToSegmentFast(p2, p3, loc));
+ double rY = sqrt(rSeg * rSeg - rX * rX);
+ // rRup is the hypoteneuse of rRup (above) and rY
+ return Distance.create(rJB, hypot(rRup, rY), rX);
+ }
+
+ double rJB = (rX > widthH) ? rX - widthH : 0.0;
+ return Distance.create(rJB, rRup, rX);
+ }
+
+ @Override
+ public double strike() {
+ return strike;
+ }
+
+ @Override
+ public double dip() {
+ return dip;
+ }
+
+ @Override
+ public double dipDirection() {
+ return (strike + 90.0) % 360.0;
+ }
+
+ @Override
+ public double length() {
+ return length;
+ }
+
+ @Override
+ public double width() {
+ return width;
+ }
+
+ @Override
+ public double area() {
+ return length * width;
+ }
+
+ @Override
+ public double depth() {
+ return zTor;
+ }
+
+ @Override
+ public Location centroid() {
+ return location;
+ }
+ }
+
+}