{-# OPTIONS -Wall #-}

{- | 
Module      :  LPFPCore.Current
Copyright   :  (c) Scott N. Walck 2023
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  stable

Code from chapter 26 of the book Learn Physics with Functional Programming
-}

module LPFPCore.Current where

import LPFPCore.SimpleVec
    ( R, Vec, sumV, (><), (*^) )
import LPFPCore.CoordinateSystems
    ( VectorField, rVF, cyl, phiHat )
import LPFPCore.Geometry
    ( Curve(..), Surface(..), Volume(..) )
import LPFPCore.ElectricField
    ( CurveApprox, curveSample, surfaceSample, volumeSample
    , vectorSurfaceIntegral, vectorVolumeIntegral )

type Current = R

data CurrentDistribution 
  = LineCurrent    Current     Curve
  | SurfaceCurrent VectorField Surface
  | VolumeCurrent  VectorField Volume
  | MultipleCurrents [CurrentDistribution]

circularCurrentLoop :: R  -- radius
                    -> R  -- current
                    -> CurrentDistribution
circularCurrentLoop :: R -> R -> CurrentDistribution
circularCurrentLoop R
radius R
i
    = R -> Curve -> CurrentDistribution
LineCurrent R
i ((R -> Position) -> R -> R -> Curve
Curve (\R
phi -> R -> R -> R -> Position
cyl R
radius R
phi R
0) R
0 (R
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi))

wireSolenoid :: R  -- radius
             -> R  -- length
             -> R  -- turns/length
             -> R  -- current
             -> CurrentDistribution
wireSolenoid :: R -> R -> R -> R -> CurrentDistribution
wireSolenoid R
radius R
len R
n R
i
    = R -> Curve -> CurrentDistribution
LineCurrent R
i ((R -> Position) -> R -> R -> Curve
Curve (\R
phi -> R -> R -> R -> Position
cyl R
radius R
phi (R
phiforall a. Fractional a => a -> a -> a
/(R
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Num a => a -> a -> a
*R
n)))
                               (-forall a. Floating a => a
piforall a. Num a => a -> a -> a
*R
nforall a. Num a => a -> a -> a
*R
len) (forall a. Floating a => a
piforall a. Num a => a -> a -> a
*R
nforall a. Num a => a -> a -> a
*R
len))

sheetSolenoid :: R  -- radius
              -> R  -- length
              -> R  -- turns/length
              -> R  -- current
              -> CurrentDistribution
sheetSolenoid :: R -> R -> R -> R -> CurrentDistribution
sheetSolenoid R
radius R
len R
n R
i
    = VectorField -> Surface -> CurrentDistribution
SurfaceCurrent (\Position
r -> (R
nforall a. Num a => a -> a -> a
*R
i) R -> Vec -> Vec
*^ VectorField
phiHat Position
r)
      (((R, R) -> Position) -> R -> R -> (R -> R) -> (R -> R) -> Surface
Surface (\(R
phi,R
z) -> R -> R -> R -> Position
cyl R
radius R
phi R
z)
       R
0 (R
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ -R
lenforall a. Fractional a => a -> a -> a
/R
2) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ R
lenforall a. Fractional a => a -> a -> a
/R
2))

wireToroid :: R  -- small radius
           -> R  -- big radius
           -> R  -- number of turns
           -> R  -- current
           -> CurrentDistribution
wireToroid :: R -> R -> R -> R -> CurrentDistribution
wireToroid R
smallR R
bigR R
n R
i
    = let alpha :: R -> R
alpha R
phi = R
n forall a. Num a => a -> a -> a
* R
phi
          curve :: R -> Position
curve R
phi = R -> R -> R -> Position
cyl (R
bigR forall a. Num a => a -> a -> a
+ R
smallR forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (R -> R
alpha R
phi)) R
phi
                      (R
smallR forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (R -> R
alpha R
phi))
      in R -> Curve -> CurrentDistribution
LineCurrent R
i ((R -> Position) -> R -> R -> Curve
Curve R -> Position
curve R
0 (R
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi))

crossedLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec
crossedLineIntegral :: CurveApprox -> VectorField -> Curve -> Vec
crossedLineIntegral CurveApprox
approx VectorField
vF Curve
c
    = [Vec] -> Vec
sumV [VectorField
vF Position
r' Vec -> Vec -> Vec
>< Vec
dl' | (Position
r',Vec
dl') <- CurveApprox
approx Curve
c]

magneticDipoleMoment :: CurrentDistribution -> Vec
magneticDipoleMoment :: CurrentDistribution -> Vec
magneticDipoleMoment (LineCurrent    R
i Curve
c)
    = CurveApprox -> VectorField -> Curve -> Vec
crossedLineIntegral   (Int -> CurveApprox
curveSample  Int
1000) (\Position
r -> R
0.5 R -> Vec -> Vec
*^ R
i R -> Vec -> Vec
*^ VectorField
rVF Position
r) Curve
c
magneticDipoleMoment (SurfaceCurrent VectorField
k Surface
s)
    = SurfaceApprox -> VectorField -> Surface -> Vec
vectorSurfaceIntegral (Int -> SurfaceApprox
surfaceSample Int
200) (\Position
r -> R
0.5 R -> Vec -> Vec
*^ (VectorField
rVF Position
r Vec -> Vec -> Vec
>< VectorField
k Position
r)) Surface
s
magneticDipoleMoment (VolumeCurrent  VectorField
j Volume
v)
    = VolumeApprox -> VectorField -> Volume -> Vec
vectorVolumeIntegral  (Int -> VolumeApprox
volumeSample   Int
50) (\Position
r -> R
0.5 R -> Vec -> Vec
*^ (VectorField
rVF Position
r Vec -> Vec -> Vec
>< VectorField
j Position
r)) Volume
v
magneticDipoleMoment (MultipleCurrents [CurrentDistribution]
ds    )
    = [Vec] -> Vec
sumV [CurrentDistribution -> Vec
magneticDipoleMoment CurrentDistribution
d | CurrentDistribution
d <- [CurrentDistribution]
ds]

helmholtzCoil :: R  -- radius
              -> R  -- current
              -> CurrentDistribution
helmholtzCoil :: R -> R -> CurrentDistribution
helmholtzCoil R
radius R
i = forall a. HasCallStack => a
undefined R
radius R
i

longStraightWire :: R  -- wire length
                 -> R  -- current
                 -> CurrentDistribution
longStraightWire :: R -> R -> CurrentDistribution
longStraightWire R
len R
i = forall a. HasCallStack => a
undefined R
len R
i

torus :: R -> R -> Surface
torus :: R -> R -> Surface
torus R
smallR R
bigR
    = ((R, R) -> Position) -> R -> R -> (R -> R) -> (R -> R) -> Surface
Surface (\(R
phi,R
alpha) -> R -> R -> R -> Position
cyl (R
bigR forall a. Num a => a -> a -> a
+ R
smallR forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos R
alpha) R
phi
                               (R
smallR forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin R
alpha))
      R
0 (R
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi) (forall a b. a -> b -> a
const R
0) (forall a b. a -> b -> a
const forall a b. (a -> b) -> a -> b
$ R
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
pi)

totalCurrent :: VectorField  -- volume current density
             -> Surface
             -> Current      -- total current through surface
totalCurrent :: VectorField -> Surface -> R
totalCurrent VectorField
j Surface
s = forall a. HasCallStack => a
undefined VectorField
j Surface
s