add SRP force

Author: jacobwilliams

README.md

@@ -27,7 +27,7 @@ The goal of the Fortran Astrodynamics Toolkit is to produce a comprehensive libr
  * Force models
    - [ ] point mass gravity field
    - [x] geopotential gravity
-   - [ ] solar radiation pressure
+   - [x] solar radiation pressure
    - [ ] atmospheric drag
    - [ ] relativistic effects
  * Reference frames
@@ -45,6 +45,7 @@ The goal of the Fortran Astrodynamics Toolkit is to produce a comprehensive libr
  * Misc
    - [x] orbital element conversions
    - [x] halo orbits
+   - [x] solar eclipsing
    - [ ] targeting and optimization
    - [ ] spacecraft engine models
 

src/lighting_module.f90

@@ -6,12 +6,12 @@
     module lighting_module
 
     use kind_module,           only: wp
-    use numbers_module,        only: pi, zero, one, two, c_light
+    use numbers_module,        only: pi, zero, one, two, c_light, fourpi, solar_luminosity
     use vector_module,         only: unit, cross, axis_angle_rotation
     use ephemeris_module,      only: ephemeris_class
     use transformation_module, only: icrf_frame
     use celestial_body_module, only: celestial_body, body_sun, body_ssb
-    use conversion_module,     only: deg2rad, rad2deg
+    use conversion_module,     only: deg2rad, rad2deg, km2m
     use math_module,           only: wrap_angle
 
     implicit none
@@ -25,12 +25,39 @@ module lighting_module
     public :: solar_fraction_alt  ! low-level routine
     public :: solar_fraction_alt2 ! low-level routine
     public :: cubic_shadow_model  ! low-level routine
+    public :: solar_radiation_pressure
 
     public :: lighting_module_test
 
     contains
 !*****************************************************************************************
 
+!********************************************************************************
+!>
+!  Compute the solar radiation pressure force vector on a spacecraft.
+
+    function solar_radiation_pressure(area, cr, r_sc_sun, sunfrac) result (srp)
+
+    real(wp),intent(in)              :: area      !! cross-sectional area of spacecraft [m^2]
+    real(wp),intent(in)              :: cr        !! coefficient of reflectivity
+    real(wp),dimension(3),intent(in) :: r_sc_sun  !! vector from spacecraft to sun [km]
+    real(wp),intent(in)              :: sunfrac   !! sun fraction [0=total eclipse, 1=no eclipse]
+    real(wp),dimension(3)            :: srp       !! solar radiation pressure force vector [N]
+
+    real(wp) :: mag   !! srp magnitude
+    real(wp) :: r_mag !! magnitude of `r_sc_sun`
+
+    if (sunfrac==zero) then
+        srp = zero
+    else
+        r_mag = norm2(r_sc_sun) * km2m ! (m)
+        mag = sunfrac * cr * solar_luminosity * area / (c_light * km2m * fourpi * r_mag**2)
+        srp = -mag * r_sc_sun * km2m / r_mag   ! (N)
+    end if
+
+    end function solar_radiation_pressure
+!********************************************************************************
+
 !********************************************************************************
 !>
 !  Compute the "sun fraction" using the selected shadow model.

src/numbers_module.f90

@@ -30,6 +30,8 @@ module numbers_module
                                                                           !! constant \( km^3/kg-s^2  \)
 
     real(wp),parameter,public :: c_light = 299792.458_wp !! speed of light in km/s
+    real(wp),parameter,public :: solar_luminosity = 3.8275e+26_wp !! solar luminosity (W)
+        !! see: "Resolution B3 on recommended nominal conversion constants for selected solar and planetary properties". IAU. 2015, which has: 3.8274 (+/- 0.0014) x 10^26 W.
 
     !> 3x3 identity matrix:
     real(wp),dimension(3,3),parameter,public :: identity_3x3 = reshape(&