diff --git a/docs/src/guide.md b/docs/src/guide.md
index ef60668..931b83d 100644
--- a/docs/src/guide.md
+++ b/docs/src/guide.md
@@ -128,6 +128,12 @@ r, c = lifting_line_geometry(grids)
 cf, cm = lifting_line_coefficients(system, r, c; frame=Body())
 nothing #hide
 ```
+Alternatively, this block of code can be simplified by calling the single overloaded `lifting_line_coefficients` function with the `system` object.
+```
+# lifting line geometry and lifting line coefficients are also combined in the single overloaded function
+cf, cm = lifting_line_coefficients(system; frame=Body())
+nothing #hide
+```
 These coefficients are defined as ``c_f = \frac{F'}{q_\infty c}`` and ``c_m = \frac{M'}{q_\infty c^2}``, respectively, where ``F'`` is the force per unit length
 along the lifting line, ``M'`` is the moment per unit length along the lifting line, ``q_\infty`` is the freestream dynamic pressure, and ``c`` is the local chord length.  
 By default, these coefficients are defined in the body frame, but may be returned in the stability or wind frame by using the `frame` keyword argument.  Note that further manipulations upon these coefficients may be required to calculate local aerodynamic coefficients since 1) the local frame of reference is not necessarily equivalent to the global frame of reference and 2) the quantities used to normalize a given local aerodynamic coefficient may vary from those used in this package.
diff --git a/src/geometry.jl b/src/geometry.jl
index d2d5837..8516b64 100644
--- a/src/geometry.jl
+++ b/src/geometry.jl
@@ -903,6 +903,96 @@ function lifting_line_geometry!(r, c, grids, xc=0.25)
     return r, c
 end
 
+"""
+    lifting_line_geometry(system, xc=0.25)
+Construct a lifting line representation of the surfaces in `system` at the
+normalized chordwise location `xc`.  Return the lifting line coordinates
+and chord lengths.
+# Arguments
+ - `system`: System containing the surfaces
+ - `xc`: Normalized chordwise location of the lifting line from the leading edge.
+    Defaults to the quarter chord
+# Return Arguments:
+ - `r`: Vector with length equal to the number of surfaces, with each element
+    being a matrix with size (3, ns+1) which contains the x, y, and z coordinates
+    of the resulting lifting line coordinates
+ - `c`: Vector with length equal to the number of surfaces, with each element
+    being a vector of length `ns+1` which contains the chord lengths at each
+    lifting line coordinate.
+"""
+function lifting_line_geometry!(system::System, xc=0.25) where TF
+    surfaces = system.surfaces
+    nsurf = length(surfaces)
+    r = Vector{Matrix{TF}}(undef, nsurf)
+    c = Vector{Vector{TF}}(undef, nsurf)
+    for isurf = 1:nsurf
+        ns = size(surfaces[isurf], 2)
+        r[isurf] = Matrix{TF}(undef, 3, ns+1)
+        c[isurf] = Vector{TF}(undef, ns+1)
+    end
+    system.lifting_line_r = r
+    system.lifting_line_c = c
+    return lifting_line_geometry!(r, c, system, xc)
+end
+
+"""
+    lifting_line_geometry!(r, c, surfaces, xc=0.25)
+In-place version of [`lifting_line_geometry`](@ref)
+"""
+function lifting_line_geometry!(r, c, system::System, xc=0.25) where TF
+    surfaces = system.surfaces
+    nsurf = length(surfaces)
+    # iterate through each lifting surface
+    for isurf = 1:nsurf
+        # extract current grid
+        surface = surfaces[isurf]
+        # dimensions of this surface
+        nc = size(surface, 1)
+        ns = size(surface, 2)
+        # loop through each spanwise section
+        for j = 1:ns
+            # get leading and trailing edges
+            le = zeros(3)
+            le .= surface[end,j].rbl
+            prev_grid = zeros(3)
+            for i = nc:-1:1
+                grid_vec = zeros(3)
+                surface_vec = zeros(3)
+                surface_vec .= surface[i,j].rtl .- surface[i,j].rbl
+                grid_vec .= (surface_vec .- prev_grid * xc) / (1 - xc)
+                prev_grid .= grid_vec
+                le .+= grid_vec
+            end
+            # le = SVector(surface[1,1,j], surface[2,1,j], surface[3,1,j])
+            te = surface[end,j].rbl
+            # get quarter-chord
+            r[isurf][:,j] = linearinterp(xc, le, te)
+            # get chord length
+            c[isurf][j] = norm(le - te)
+        end
+
+        le = zeros(3)
+        le .= surface[end,end].rbr
+        prev_grid = zeros(3)
+        for i = nc:-1:1
+            grid_vec = zeros(3)
+            surface_vec = zeros(3)
+            surface_vec .= surface[i,end].rtr .- surface[i,end].rbr
+            grid_vec .= (surface_vec .- prev_grid * xc) / (1 - xc)
+            prev_grid .= grid_vec
+            le .+= grid_vec
+        end
+        # le = SVector(surface[1,1,end], surface[2,1,end], surface[3,1,end])
+        te = surface[end,end].rbr
+        # get quarter-chord
+        r[isurf][:,end] = linearinterp(xc, le, te)
+        # get chord length
+        c[isurf][end] = norm(le - te)
+        
+    end
+    return r, c
+end
+
 """
     translate(grid, r)
 
@@ -1015,4 +1105,4 @@ end
 
 Test whether and of the points in `args` are not on the symmetry plane (y = 0)
 """
-not_on_symmetry_plane(args...; tol=eps()) = !on_symmetry_plane(args...; tol=tol)
+not_on_symmetry_plane(args...; tol=eps()) = !on_symmetry_plane(args...; tol=tol)
\ No newline at end of file
diff --git a/src/nearfield.jl b/src/nearfield.jl
index ffcae90..0e4edb0 100644
--- a/src/nearfield.jl
+++ b/src/nearfield.jl
@@ -914,6 +914,45 @@ function lifting_line_coefficients(system, r, c; frame=Body())
     return lifting_line_coefficients!(cf, cm, system, r, c; frame)
 end
 
+"""
+    lifting_line_coefficients(system; frame=Body())
+
+Return the force and moment coefficients (per unit span) for each spanwise segment
+of a lifting line representation of the geometry.
+
+This function requires that a near-field analysis has been performed on `system`
+to obtain panel forces.
+
+# Arguments
+ - `system`: Object of type [`System`](@ref) that holds precalculated
+    system properties.
+
+# Keyword Arguments
+ - `frame`: frame in which to return `cf` and `cm`, possible options are
+    [`Body()`](@ref) (default), [`Stability()`](@ref), and [`Wind()`](@ref)`
+
+# Return Arguments:
+ - `cf`: Vector with length equal to the number of surfaces, with each element
+    being a matrix with size (3, ns) which contains the x, y, and z direction
+    force coefficients (per unit span) for each spanwise segment.
+ - `cm`: Vector with length equal to the number of surfaces, with each element
+    being a matrix with size (3, ns) which contains the x, y, and z direction
+    moment coefficients (per unit span) for each spanwise segment.
+"""
+function lifting_line_coefficients(system; frame=Body())
+    lifting_line_geometry!(system)
+    TF = promote_type(eltype(system), eltype(eltype(r)), eltype(eltype(c)))
+    nsurf = length(system.surfaces)
+    cf = Vector{Matrix{TF}}(undef, nsurf)
+    cm = Vector{Matrix{TF}}(undef, nsurf)
+    for isurf = 1:nsurf
+        ns = size(system.surfaces[isurf], 2)
+        cf[isurf] = Matrix{TF}(undef, 3, ns)
+        cm[isurf] = Matrix{TF}(undef, 3, ns)
+    end
+    return lifting_line_coefficients!(cf, cm, system; frame)
+end
+
 """
     lifting_line_coefficients!(cf, cm, system, r, c; frame=Body())
 
@@ -981,6 +1020,75 @@ function lifting_line_coefficients!(cf, cm, system, r, c; frame=Body())
     return cf, cm
 end
 
+"""
+    lifting_line_coefficients!(cf, cm, system; frame=Body())
+
+In-place version of [`lifting_line_coefficients`](@ref)
+"""
+function lifting_line_coefficients!(cf, cm, system; frame=Body())
+    r = system.lifting_line_r
+    c = system.lifting_line_c
+
+    # number of surfaces
+    nsurf = length(system.surfaces)
+
+    # check that a near field analysis has been performed
+    @assert system.near_field_analysis[] "Near field analysis required"
+
+    # extract reference properties
+    ref = system.reference[]
+    fs = system.freestream[]
+
+    # iterate through each lifting surface
+    for isurf = 1:nsurf
+        nc, ns = size(system.surfaces[isurf])
+        # extract current surface panels and panel properties
+        panels = system.surfaces[isurf]
+        properties = system.properties[isurf]
+        # loop through each chordwise set of panels
+        for j = 1:ns
+            # calculate segment length
+            rls = SVector(r[isurf][1,j], r[isurf][2,j], r[isurf][3,j])
+            rrs = SVector(r[isurf][1,j+1], r[isurf][2,j+1], r[isurf][3,j+1])
+            ds = norm(rrs - rls)
+            # calculate reference location
+            rs = (rls + rrs)/2
+            # calculate reference chord
+            cs = (c[isurf][j] + c[isurf][j+1])/2
+            # calculate section force and moment coefficients
+            cfj = @SVector zeros(eltype(cf[isurf]), 3)
+            cmj = @SVector zeros(eltype(cm[isurf]), 3)
+            for i = 1:nc
+                # add influence of bound vortex
+                rb = top_center(panels[i,j])
+                cfb = properties[i,j].cfb
+                cfj += cfb
+                cmj += cross(rb - rs, cfb)
+                # add influence of left vortex leg
+                rl = left_center(panels[i,j])
+                cfl = properties[i,j].cfl
+                cfj += cfl
+                cmj += cross(rl - rs, cfl)
+                # add influence of right vortex leg
+                rr = right_center(panels[i,j])
+                cfr = properties[i,j].cfr
+                cfj += cfr
+                cmj += cross(rr - rs, cfr)
+            end
+            # update normalization
+            cfj *= ref.S/(ds*cs)
+            cmj *= ref.S/(ds*cs^2)
+            # change coordinate frame
+            cfj, cmj = body_to_frame(cfj, cmj, ref, fs, frame)
+            # save coefficients
+            cf[isurf][:,j] = cfj
+            cm[isurf][:,j] = cmj
+        end
+    end
+
+    return cf, cm
+end
+
 """
     body_to_frame(CF, CM, reference, freestream, frame)
 
diff --git a/src/system.jl b/src/system.jl
index ad06a3f..c48bf5e 100644
--- a/src/system.jl
+++ b/src/system.jl
@@ -103,6 +103,8 @@ struct System{TF}
     Vv::Vector{Matrix{SVector{3, TF}}}
     Vte::Vector{Vector{SVector{3, TF}}}
     dΓdt::Vector{TF}
+    lifting_line_r::Vector{Matrix{TF}}
+    lifting_line_c::Vector{Vector{TF}}
 end
 
 @inline Base.eltype(::Type{System{TF}}) where TF = TF
diff --git a/test/runtests.jl b/test/runtests.jl
index 48c492b..631b410 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -908,21 +908,27 @@ end
 
     system = steady_analysis(surfaces, ref, fs; symmetric=symmetric)
 
+    # old method to generate lifting line coefficients
     r_ll, c_ll = lifting_line_geometry(grids)
-
     cf, cm = lifting_line_coefficients(system, r_ll, c_ll; frame=Stability())
 
+    # new overloaded function that calls lifting_line_geometry
+    alt_cf, alt_cm = lifting_line_coefficients(system; frame=Stability())
+
     cl_avl = [0.2618, 0.2646, 0.2661, 0.2664, 0.2654, 0.2628, 0.2584, 0.2513,
         0.2404, 0.2233, 0.1952, 0.1434]
     cd_avl = [0.0029, 0.0024, 0.0023, 0.0023, 0.0023, 0.0023, 0.0024, 0.0024,
         0.0025, 0.0026, 0.0026, 0.0022]
     cm_avl = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 
+    @test isapprox(cf, alt_cf, atol=1e-3, norm=(x)->norm(x, Inf))
+    @test isapprox(cm, alt_cm, atol=1e-3, norm=(x)->norm(x, Inf))
     @test isapprox(cf[1][3,:], cl_avl, atol=1e-3, norm=(x)->norm(x, Inf))
     @test isapprox(cf[1][1,:], cd_avl, atol=1e-4, norm=(x)->norm(x, Inf))
     @test isapprox(cm[1][2,:], cm_avl, atol=1e-4, norm=(x)->norm(x, Inf))
 end
 
+
 @testset "Geometry Generation" begin
 
     # Tests of the geometry generation functions