Pratt and Chandler’s theory calculation

gistfile1.hs

import Data.List
import Data.Maybe
import Control.Monad

type RealFunc = Double -> Double

fi :: Int -> Double
fi = fromIntegral

integral :: RealFunc -> (Double,Double) -> Maybe Int -> Double
integral f (from,to) m_num_interval = sum (map f xs) * dx
  where
    xs = map (\i -> from + (fi i+0.5)*dx) [0..num_interval-1]
    dx = (to-from) / fi num_interval
    num_interval = fromMaybe 10000 m_num_interval

func_j :: Double -> Int -> Int -> RealFunc
func_j r m n x = ((x-r)/r)** fi (m+n) * (4*pi*x**2)

int_j :: Env -> Int -> Int -> Double
int_j env m n
  | m+n == 0 = 4*pi*(r_max env)**3 / 3
  | otherwise = integral (func_j (r_max env) m n) ((r_min env),r_max env) Nothing

phi :: Int -> Double -> RealFunc
phi n r x | x >= r = 0
        | otherwise = ((x-r)/r)**fi n

phi_hat_0,phi_hat_1,phi_hat_2,phi_hat_3 :: Double -> RealFunc

phi_hat_0 r k = pi*4.0 * (sin kr / k**3 - r * cos kr / (k**2))
  where kr = k * r
phi_hat_1 r k = pi*4.0 * (sin kr / k**3 + 2.0 * (cos kr-1)/ (r*k**4))
  where kr = k * r

phi_hat_2 r k = pi*4.0 * (-6.0) *(sin kr / (r**2*k**5) + (2 * r * cos kr + 4) / (r*k**4))
  where kr = k * r

phi_hat_3 r k = pi*4.0 * (6*sin kr / (r**2*k**5) + 24 * (1-cos kr)/(r**3*k**6) - 6/(r*k**4))
  where kr = k * r

phi_hat :: Int -> (Double -> RealFunc)
phi_hat 0 = phi_hat_0
phi_hat 1 = phi_hat_1
phi_hat 2 = phi_hat_2
phi_hat 3 = phi_hat_3

func_k :: Double -> Int -> Int -> RealFunc
func_k r m n k = k**2 * (phi_hat m r k) * (phi_hat n r k) * s_minus_1 k

s_minus_1 :: RealFunc
s_minus_1 x = ya + (x-xa)*(yb-ya)/(xb-xa)
  where
    (as,bs) = span ((<= x) . fst) s_minus_1_table
    ((xa,ya),(xb,yb)) = (last as,head bs)


int_k :: Env -> Int -> Int -> Double
int_k env m n = (4*pi/(2*pi)**3) * integral (func_k (r_max env) m n) (r_min env, r_max env) Nothing

data Env = Env{
  r_min :: Double,
  r_max :: Double
}

defEnv = Env 0 2.7

mat_a env = map j_plus_k [(x,y) | x <- [0..3], y <- [0..3]]
  where j_plus_k (m,n) = ((m,n),int_j env m n + int_k env m n)

int_i :: Env -> Int -> Double
int_i env n = integral (phi n (r_max env)) ((r_min env),(r_max env)) Nothing

func_i r n x = 4*pi*x**2*(phi r n x)

vec_i env = map (int_i env) [0..3]

print_mat m n mat = do
  putStr "["
  forM_ [0..m-1] $ \i -> do
    putStr $ intercalate "\t" $ map show $ map (maybe (0/0) snd) $ map (\j -> find ( (== (i,j)) . fst) mat) [0..n-1]
    putStrLn ";"
  putStrLn "]"

main = do
  let envs = map (Env 0.1) $ map (\i -> 3.0+fi i*1.0) [0..2]
  forM_ [0..2] $ \i -> do
    let env = envs !! i
    putStrLn $ "r_max = " ++ show (r_max env)
    putStr $ "A{" ++ show (i+1) ++ "} = "
    print_mat 4 4 (mat_a env)
    putStr $ "I{" ++ show (i+1) ++ "} = "
    print $ vec_i env

s_minus_1_table :: [(Double,Double)]
s_minus_1_table =
  [( 0.0,-0.938),
      ( 0.08,-0.938),
      ( 0.16,-0.938),
      ( 0.24,-0.937),
      ( 0.32,-0.937),
      ( 0.4,-0.936),
      ( 0.48,-0.935),
      ( 0.56,-0.932),
      ( 0.64,-0.928),
      ( 0.72,-0.921),
      ( 0.8,-0.913),
      ( 0.88,-0.902),
      ( 0.96,-0.890),
      ( 1.04,-0.878),
      ( 1.12,-0.861),
      ( 1.2,-0.835),
      ( 1.28,-0.796),
      ( 1.36,-0.748),
      ( 1.44,-0.685),
      ( 1.52,-0.599),
      ( 1.6,-0.478),
      ( 1.68,-0.326),
      ( 1.76,-0.158),
      ( 1.84,-0.002),
      ( 1.92,0.119),
      ( 2.00,0.190),
      ( 2.08,0.212),
      ( 2.16,0.197),
      ( 2.24,0.166),
      ( 2.32,0.134),
      ( 2.4,0.116),
      ( 2.48,0.119),
      ( 2.56,0.147),
      ( 2.64,0.193),
      ( 2.72,0.245),
      ( 2.8,0.287),
      ( 2.88,0.301),
      ( 2.96,0.277),
      ( 3.04,0.215),
      ( 3.12,0.125),
      ( 3.20,0.022),
      ( 3.28,-0.079),
      ( 3.36,-0.164),
      ( 3.44,-0.227),
      ( 3.52,-0.262),
      ( 3.6,-0.271),
      ( 3.68,-0.255),
      ( 3.76,-0.233),
      ( 3.84,-0.199),
      ( 3.92,-0.162),
      ( 4.0,-0.124),
      ( 4.08,-0.083),
      ( 4.16,-0.039),
      ( 4.24,0.005),
      ( 4.32,0.045),
      ( 4.4,0.076),
      ( 4.48,0.098),
      ( 4.56,0.114),
      ( 4.64,0.126),
      ( 4.72,0.135),
      ( 4.8,0.14),
      ( 4.88,0.136),
      ( 4.96,0.124),
      ( 5.04,0.104  ),
      ( 5.12,0.080),
      ( 5.2,0.057),
      ( 5.28,0.035),
      ( 5.36,0.013),
      ( 5.44,-0.009),
      ( 5.52,-0.034),
      ( 5.6,-0.057),
      ( 5.68,-0.078),
      ( 5.76,-0.093),
      ( 5.84,-0.102),
      ( 5.92,-0.105),
      ( 6.0,-0.102),
      ( 6.08,-0.094),
      ( 6.16,-0.080),
      ( 6.24,-0.063),
      ( 6.32,-0.044),
      ( 6.4,-0.026),
      ( 6.48,-0.009),
      ( 6.56,0.007),
      ( 6.64,0.023),
      ( 6.72,0.040),
      ( 6.8,0.056),
      ( 6.88,0.069),
      ( 6.96,0.075),
      ( 7.04,0.077),
      ( 7.12,0.074),
      ( 7.2,0.069),
      ( 7.28,0.063),
      ( 7.36,0.053),
      ( 7.44,0.040),
      ( 7.52,0.024),
      ( 7.6,0.007),
      ( 7.68,-0.008),
      ( 7.76,-0.020),
      ( 7.84,-0.030),
      ( 7.92,-0.037),
      ( 8.0,-0.044),
      ( 8.08,-0.050),
      ( 8.16,-0.054),
      ( 8.24,-0.054),
      ( 8.32,-0.050),
      ( 8.4,-0.042),
      ( 8.48,-0.033),
      ( 8.56,-0.023),
      ( 8.64,-0.013),
      ( 8.72,-0.005),
      ( 8.8,0.003),
      ( 8.88,0.010),
      ( 8.96,0.017),
      ( 9.04,0.023),
      ( 9.12,0.028),
      ( 9.2,0.031),
      ( 9.28,0.032),
      ( 9.36,0.031),
      ( 9.44,0.029),
      ( 9.52,0.027),
      ( 9.6,0.025),
      ( 9.68,0.021),
      ( 9.76,0.015),
      ( 9.84,0.007),
      ( 9.92,-0.001),
      ( 10.00,-0.009),
      ( 10.08,-0.014),
      ( 10.16,-0.017),
      ( 10.24,-0.018),
      ( 10.32,-0.019),
      ( 10.4,-0.020),
      ( 10.48,-0.020),
      ( 10.56,-0.020),
      ( 10.64,-0.019),
      ( 10.72,-0.016),
      ( 10.8,-0.012),
      ( 10.88,-0.009),
      ( 10.96,-0.007),
      ( 11.04,-0.005),
      ( 11.12,-0.003),
      ( 11.20,-0.0),
      ( 11.28,0.003),
      ( 11.36,0.007),
      ( 11.44,0.010),
      ( 11.52,0.013),
      ( 11.6,0.016),
      ( 11.68,0.017),
      ( 11.76,0.017),
      ( 11.84,0.015),
      ( 11.92,0.013),
      ( 12.0,0.010),
      ( 12.08,0.008),
      ( 12.16,0.006),
      ( 12.24,0.003),
      ( 12.32,0.000),
      ( 12.4,-0.002),
      ( 12.48,-0.004),
      ( 12.56,-0.005),
      ( 12.64,-0.005),
      ( 12.72,-0.006),
      ( 12.8,-0.007),
      ( 12.88,-0.010),
      ( 12.96,-0.012),
      ( 13.04,-0.011),
      ( 13.12,-0.009),
      ( 13.2,-0.006),
      ( 13.28,-0.003),
      ( 13.36,-0.001),
      ( 13.44,0.001),
      ( 13.52,0.002),
      ( 13.6,0.005),
      ( 13.68,0.007),
      ( 13.76,0.008),
      ( 13.84,0.008),
      ( 13.92,0.007),
      ( 14.0,0.005),
      ( 14.08,0.004),
      ( 14.16,0.004),
      ( 14.24,0.004),
      ( 14.32,0.004),
      ( 14.4,0.002),
      ( 14.48,0.000),
      ( 14.56,-0.002),
      ( 14.64,-0.003),
      ( 14.72,-0.004),
      ( 14.8,-0.004),
      ( 14.88,-0.004),
      ( 14.96,-0.004),
      ( 15.04,-0.004),
      ( 15.12,-0.003),
      ( 15.2,-0.003),
      ( 15.28,-0.003),
      ( 15.36,-0.003),
      ( 15.44,-0.003),
      ( 15.52,-0.002),
      ( 15.6,-0.001),
      ( 15.68,-0.000),
      ( 15.76,-0.000),
      ( 15.84,-0.001),
      ( 15.92,-0.002),
      ( 16.0,-0.003)]