using MieSeries
using Base.Test


p = [1.0, 2.1, -3.7]
K = 4

Y = MieSeries.Ylm(K,p)
Y_from_matlab = [
    0.282094791773878 + 0.000000000000000im
    0.079054029942852 - 0.166013462879989im
    -0.413657340830090 + 0.000000000000000im
    -0.079054029942852 - 0.166013462879989im
    -0.068963090518276 - 0.084939876884680im
    -0.149655973558723 + 0.314277544473318im
    0.362782863277375 + 0.000000000000000im
    0.149655973558723 + 0.314277544473318im
    -0.068963090518276 + 0.084939876884680im
    -0.061128748571790 + 0.014799854825926im
    0.154472358060769 + 0.190259209341709im
    0.191065304223659 - 0.401237138869684im
    -0.184433850338810 + 0.000000000000000im
    -0.191065304223659 - 0.401237138869684im
    0.154472358060769 - 0.190259209341709im
    0.061128748571790 + 0.014799854825926im]

norm(Y-Y_from_matlab)
for i in eachindex(Y)
    @eval @test_approx_eq Y[$i] Y_from_matlab[$i]
end

Xx, Xy, Xz, Yx, Yy, Yz, Zx, Zy, Zz = MieSeries.Xlm(K,p)


Xx_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    -0.206828670415045 + 0.000000000000000im
    0.000000000000000 - 0.166013462879989im
    -0.206828670415045 + 0.000000000000000im
    -0.061096795363054 + 0.128303270262413im
    0.153237367829497 - 0.034676559530382im
    0.000000000000000 + 0.314277544473318im
    0.153237367829497 + 0.034676559530382im
    0.061096795363054 + 0.128303270262413im
    0.054614225945323 + 0.067266788554357im
    0.065596704579624 - 0.177906321035805im
    -0.021710262651990 + 0.086841050607960im
    0.000000000000000 - 0.401237138869684im
    -0.021710262651990 - 0.086841050607960im
    -0.065596704579624 - 0.177906321035805im
    0.054614225945323 - 0.067266788554357im]

Xy_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    0.000000000000000 + 0.206828670415045im
    0.000000000000000 + 0.079054029942852im
    0.000000000000000 - 0.206828670415045im
    0.128303270262413 + 0.061096795363054im
    0.034676559530382 - 0.209545495447878im
    0.000000000000000 - 0.149655973558723im
    0.034676559530382 + 0.209545495447878im
    -0.128303270262413 + 0.061096795363054im
    0.067266788554357 - 0.054614225945323im
    -0.188371398743793 - 0.108821257220184im
    -0.086841050607960 + 0.162723587686820im
    0.000000000000000 + 0.191065304223659im
    -0.086841050607960 - 0.162723587686820im
    0.188371398743793 - 0.108821257220184im
    0.067266788554357 + 0.054614225945323im]

Xz_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    -0.055899640652715 + 0.117389245370701im
    0.000000000000000 + 0.000000000000000im
    -0.055899640652715 - 0.117389245370701im
    0.056308127618382 + 0.069353119060764im
    0.061096795363054 - 0.128303270262413im
    0.000000000000000 + 0.000000000000000im
    0.061096795363054 + 0.128303270262413im
    -0.056308127618382 + 0.069353119060764im
    0.052939049164722 - 0.012817050251573im
    -0.089184657508741 - 0.109846205729241im
    -0.055155802413164 + 0.115827185067644im
    0.000000000000000 + 0.000000000000000im
    -0.055155802413164 - 0.115827185067644im
    0.089184657508741 - 0.109846205729241im
    0.052939049164722 + 0.012817050251573im]

@test norm(Xx - Xx_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Xy - Xy_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Xz - Xz_from_matlab) < sqrt(eps(eltype(p)))


Yx_from_matlab = [
    0.064547345156897 + 0.000000000000000im
    0.018088698925201 - 0.037986267742922im
    -0.094650748379100 + 0.000000000000000im
    -0.018088698925201 - 0.037986267742922im
    -0.015779746867279 - 0.019435465349727im
    -0.034243438949520 + 0.071911221793991im
    0.083009936290289 + 0.000000000000000im
    0.034243438949520 + 0.071911221793991im
    -0.015779746867279 + 0.019435465349727im
    -0.013987136764422 + 0.003386419620560im
    0.035345496987908 + 0.043534043210913im
    0.043718489312336 - 0.091808827555906im
    -0.042201117296689 + 0.000000000000000im
    -0.043718489312336 - 0.091808827555906im
    0.035345496987908 - 0.043534043210913im
    0.013987136764422 + 0.003386419620560im]

Yy_from_matlab = [
    0.135549424829483 + 0.000000000000000im
    0.037986267742922 - 0.079771162260136im
    -0.198766571596109 + 0.000000000000000im
    -0.037986267742922 - 0.079771162260136im
    -0.033137468421285 - 0.040814477234427im
    -0.071911221793991 + 0.151013565767382im
    0.174320866209607 + 0.000000000000000im
    0.071911221793991 + 0.151013565767382im
    -0.033137468421285 + 0.040814477234427im
    -0.029372987205287 + 0.007111481203177im
    0.074225543674607 + 0.091421490742917im
    0.091808827555906 - 0.192798537867402im
    -0.088622346323047 + 0.000000000000000im
    -0.091808827555906 - 0.192798537867402im
    0.074225543674607 - 0.091421490742917im
    0.029372987205287 + 0.007111481203177im]

Yz_from_matlab = [
    -0.238825177080517 + 0.000000000000000im
    -0.066928186023243 + 0.140549190648810im
     0.350207769002669 + 0.000000000000000im
     0.066928186023243 + 0.140549190648810im
     0.058385063408931 + 0.071911221793991im
     0.126700724113223 - 0.266071520637768im
    -0.307136764274069 + 0.000000000000000im
    -0.126700724113223 - 0.266071520637768im
     0.058385063408931 - 0.071911221793991im
     0.051752406028363 - 0.012529752596074im
    -0.130778338855260 - 0.161075959880378im
    -0.161758410455643 + 0.339692661956851im
     0.156144133997750 + 0.000000000000000im
     0.161758410455643 + 0.339692661956851im
    -0.130778338855260 + 0.161075959880378im
    -0.051752406028363 - 0.012529752596074im]

@test norm(Yx - Yx_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Yy - Yy_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Yz - Yz_from_matlab) < sqrt(eps(eltype(p)))

Zx_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    -0.026860347512988 + 0.231510614278609im
    0.000000000000000 + 0.066928186023243im
    -0.026860347512988 - 0.231510614278609im
    0.135679873108069 + 0.085050235134019im
    0.058715266725129 - 0.239055014523741im
    0.000000000000000 - 0.126700724113223im
    0.058715266725129 + 0.239055014523741im
    -0.135679873108069 + 0.085050235134019im
    0.082386704863588 - 0.052395848498920im
    -0.202331852112205 - 0.144911736140107im
    -0.100023706233593 + 0.193420082720831im
    0.000000000000000 + 0.161758410455644im
    -0.100023706233593 - 0.193420082720831im
    0.202331852112205 - 0.144911736140107im
    0.082386704863588 + 0.052395848498920im]

Zy_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    0.187894526174186 - 0.026860347512988im
    0.000000000000000 + 0.140549190648810im
    0.187894526174186 + 0.026860347512988im
    0.038841244654719 - 0.124492234448281im
    -0.143712605222459 + 0.058715266725129im
    0.000000000000000 - 0.266071520637768im
    -0.143712605222459 - 0.058715266725129im
    -0.038841244654719 - 0.124492234448281im
    -0.058350340882630 - 0.054016226339494im
    -0.035128304403982 + 0.175752233030384im
    0.031000635931996 - 0.100023706233593im
    0.000000000000000 + 0.339692661956851im
    0.031000635931996 + 0.100023706233593im
    0.035128304403982 + 0.175752233030384im
    -0.058350340882630 + 0.054016226339494im]

Zz_from_matlab = [
    0.000000000000000 + 0.000000000000000im
    0.099383285798055 + 0.047325374189550im
    0.000000000000000 + 0.097859861185336im
    0.099383285798055 - 0.047325374189550im
    0.058715266725129 - 0.047671204650641im
    -0.065697622768118 - 0.031284582270532im
    0.000000000000000 - 0.185257004716901im
    -0.065697622768118 + 0.031284582270532im
    -0.058715266725129 - 0.047671204650641im
    -0.010851084051334 - 0.044818898327529im
    -0.074621970637991 + 0.060585933303702im
    -0.009438478588217 - 0.004494513613437im
    0.000000000000000 + 0.236517027179738im
    -0.009438478588217 + 0.004494513613437im
    0.074621970637991 + 0.060585933303702im
    -0.010851084051334 + 0.044818898327529im]

@test norm(Zx - Zx_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Zy - Zy_from_matlab) < sqrt(eps(eltype(p)))
@test norm(Zz - Zz_from_matlab) < sqrt(eps(eltype(p)))
