薄膜干渉のｼﾐｭﾚｰｼｮﾝするやつ

plot.sh

#!/bin/sh

rm -rf ./data
mkdir -p ./data

runghc ./two-beam-interference.hs
runghc ./TBI-with-lambda.hs
runghc ./TBI-with-d.hs
runghc ./TBI-with-n.hs

deg=`ls -x data/*degree.data | sed 's/  /\",\"/g' | sed 's/^/\"/g' | sed 's/$/\"/g' | tr -d "\n" | sed 's/\"\"/\",\"/g'`

lamb=`ls -x data/*nm.data | sed 's/  /\",\"/g' | sed 's/^/\"/g' | sed 's/$/\"/g' | tr -d "\n" | sed 's/\"\"/\",\"/g'`

d=`ls -x data/d*microm.data | sed 's/  /\",\"/g' | sed 's/^/\"/g' | sed 's/$/\"/g' | tr -d "\n" | sed 's/\"\"/\",\"/g'`

n=`ls -x data/n*.data | sed 's/[[:blank:]]/\",\"/g' | sed 's/^/\"/g' | sed 's/$/\"/g' | tr -d "\n" | sed 's/\"\"/\",\"/g'`

gnuplot <<EOF
    set title "Two-beam Interference"
    set xlabel "wave length[nm]"
    set ylabel "reflectance"
    set key opaque box
    set key right bottom
    set terminal png
    set output "output_deg.png"
    plot ${deg}
    set terminal pdf
    set output "output_deg.pdf"
    replot
EOF

gnuplot <<EOF
    set title "Two-beam Interference"
    set xlabel "degrees"
    set ylabel "reflectance"
    set key opaque box
    set key right bottom
    set terminal png
    set output "output_nm.png"
    plot ${lamb}
    set terminal pdf
    set output "output_nm.pdf"
    replot
EOF

gnuplot <<EOF
    set title "Two-beam Interference"
    set xlabel "wave length[nm]"
    set ylabel "reflectance"
    set key opaque box
    set key right bottom 
    set terminal png
    set output "output_delta.png"
    plot ${d}
    set terminal pdf
    set output "output_delta.pdf"
    replot
EOF

gnuplot <<EOF
    set title "Two-beam Interference"
    set xlabel "wave length[nm]"
    set ylabel "reflectance"
    set key opaque box
    set key right bottom 
    set terminal png
    set output "output_n.pdf"
    plot ${n}
    set terminal pdf
    set output "output_n.pdf"
    replot
EOF

TBI-with-d.hs

{-
波長lambdaは380nm-1500nm程度として計算する。
媒質の屈折率n1は1.5として考える。
-}

import System.IO

main :: IO ()
main = mapM_ simulate [0.1,0.3,0.4,0.5,0.6,0.7,1.0]

simulate :: Float -> IO ()
simulate ind = do
        outh <- openFile ("data/d" ++ show ind ++ "microm.data") WriteMode
        mapM_ (hPutStrLn outh) $ map (irradiance izero r d thetai) [380..780]
        hClose outh
    where
        ni = 1.0
        nt = 1.5
        r = reflectionS ni nt thetai thetat
        thetai' deg = deg / 180 * pi
        thetai = thetai' 0
        thetat = asin $ ni * sin thetai / nt
        d = ind * 0.000001
        izero = 1

irradiance :: Float -> Float -> Float -> Float -> Float -> [Char]
irradiance izero r d theta l = show l ++ " " ++ show irrad
    where
        rt = 1 - r^2
        beta = 2 * pi * d * cos theta / l * 10^9
        irrad = izero * r^2 * ( 1 + rt^2 - 2 * rt * cos (2 * beta) )

reflectionP :: Float -> Float -> Float -> Float -> Float 
reflectionP ni nt thetai thetat = (nt * cos thetai - ni * cos thetat) / (ni * cos thetat + nt * cos thetai)

reflectionS :: Float -> Float -> Float -> Float -> Float 
reflectionS ni nt thetai thetat = (ni * cos thetai - nt * cos thetat) / (ni * cos thetat + nt * cos thetai)

TBI-with-lambda.hs

{-
Two Beam Interference with lambda
波長lambdaは380nm-780nm程度として計算する。
媒質の屈折率n1は1.5として考える。
-}

import System.IO

main :: IO ()
main = mapM_ simulate [420,500,670]

simulate :: Float -> IO ()
simulate inlamb = do
        outh <- openFile ("data/" ++ show inlamb ++ "nm.data") WriteMode
        mapM_ (hPutStrLn outh) $ map (irradiance izero d inlamb) [0..70]
        hClose outh
    where
        d = 0.000001
        izero = 1

irradiance :: Float -> Float -> Float -> Float -> [Char]
irradiance izero d l theta = show theta ++ " " ++ show irrad
    where
        ni = 1.0
        nt = 1.5
        thetai = theta / 180 * pi
        r = reflectionS ni nt thetai thetat
        thetat = asin $ ni * sin thetai / nt
        rt = 1 - r^2
        beta = 2 * pi * d * cos thetai / l * 10^9
        irrad = izero * r^2 * ( 1 + rt^2 - 2 * rt * cos (2 * beta) )

reflectionP :: Float -> Float -> Float -> Float -> Float 
reflectionP ni nt thetai thetat = (nt * cos thetai - ni * cos thetat) / (ni * cos thetat + nt * cos thetai)

reflectionS :: Float -> Float -> Float -> Float -> Float 
reflectionS ni nt thetai thetat = (ni * cos thetai - nt * cos thetat) / (ni * cos thetat + nt * cos thetai)

TBI-with-n.hs

{-
波長lambdaは380nm-1500nm程度として計算する。
-}

import System.IO

main :: IO ()
main = mapM_ simulate [1.5,2.0,2.5,3.0,4.0,5.0,8.0]

simulate :: Float -> IO ()
simulate inn = do
        outh <- openFile ("data/n" ++ show inn ++ ".data") WriteMode
        mapM_ (hPutStrLn outh) $ map (irradiance izero r d thetai) [380..780]
        hClose outh
    where
        ni = 1.0
        nt = inn
        r = reflectionS ni nt thetai thetat
        thetai' deg = deg / 180 * pi
        thetai = thetai' 0
        thetat = asin $ ni * sin thetai / nt
        d = 0.000001
        izero = 1

irradiance :: Float -> Float -> Float -> Float -> Float -> [Char]
irradiance izero r d theta l = show l ++ " " ++ show irrad
    where
        rt = 1 - r^2
        beta = 2 * pi * d * cos theta / l * 10^9
        irrad = izero * r^2 * ( 1 + rt^2 - 2 * rt * cos (2 * beta) )

reflectionP :: Float -> Float -> Float -> Float -> Float 
reflectionP ni nt thetai thetat = (nt * cos thetai - ni * cos thetat) / (ni * cos thetat + nt * cos thetai)

reflectionS :: Float -> Float -> Float -> Float -> Float 
reflectionS ni nt thetai thetat = (ni * cos thetai - nt * cos thetat) / (ni * cos thetat + nt * cos thetai)

two-beam-interference.hs

{-
波長lambdaは380nm-1500nm程度として計算する。
媒質の屈折率n1は1.5として考える。
-}

import System.IO

main :: IO ()
main = mapM_ simulate [10,15..70]

simulate :: Float -> IO ()
simulate indeg = do
        outh <- openFile ("data/" ++ show indeg ++ "degree.data") WriteMode
        mapM_ (hPutStrLn outh) $ map (irradiance izero r d thetai) [380..780]
        hClose outh
    where
        ni = 1.0
        nt = 1.5
        r = reflectionS ni nt thetai thetat
        thetai' deg = deg / 180 * pi
        thetai = thetai' indeg
        thetat = asin $ ni * sin thetai / nt
        d = 0.000001
        izero = 1

irradiance :: Float -> Float -> Float -> Float -> Float -> [Char]
irradiance izero r d theta l = show l ++ " " ++ show irrad
    where
        rt = 1 - r^2
        beta = 2 * pi * d * cos theta / l * 10^9
        irrad = izero * r^2 * ( 1 + rt^2 - 2 * rt * cos (2 * beta) )

reflectionP :: Float -> Float -> Float -> Float -> Float 
reflectionP ni nt thetai thetat = (nt * cos thetai - ni * cos thetat) / (ni * cos thetat + nt * cos thetai)

reflectionS :: Float -> Float -> Float -> Float -> Float 
reflectionS ni nt thetai thetat = (ni * cos thetai - nt * cos thetat) / (ni * cos thetat + nt * cos thetai)