# File   : spherical_cylinder_approx.jl
# License: GNU v3.0
# Author : Andrei Leonard Nicusan <a.l.nicusan@bham.ac.uk>
# Date   : 28.08.2021


using DelimitedFiles
using LinearAlgebra

using GLMakie
using GeometryBasics


# Vertices and face indices generated using AutoDesk; two spherical caps were created, then the
# convex hull was saved in the "spherical_cylinder.obj" file
verts = [
    -0.945817 0.324700 0.000000
    -0.947808 0.078538 -0.309017
    -1.000000 -0.000000 0.000000
    -0.921954 -0.233471 -0.309017
    -0.945817 -0.324700 0.000000
    -0.796193 -0.520179 -0.309017
    -0.789141 -0.614213 0.000000
    -0.584151 -0.750517 -0.309017
    -0.546948 -0.837166 -0.000000
    -0.308807 -0.899526 -0.309017
    -0.245485 -0.969400 -0.000000
    0.000000 -0.951056 -0.309017
    0.082579 -0.996584 -0.000000
    0.308807 -0.899526 -0.309017
    0.401695 -0.915773 -0.000000
    0.584151 -0.750517 -0.309017
    0.677282 -0.735724 -0.000000
    0.796193 -0.520179 -0.309017
    0.879474 -0.475947 0.000000
    0.921954 -0.233471 -0.309017
    0.986361 -0.164595 0.000000
    0.947808 0.078538 -0.309017
    0.986361 0.164595 0.000000
    0.870952 0.382035 -0.309017
    0.879474 0.475947 0.000000
    0.699715 0.644133 -0.309017
    0.677282 0.735724 0.000000
    0.452653 0.836429 -0.309017
    0.401695 0.915773 0.000000
    0.156539 0.938085 -0.309017
    0.082579 0.996584 0.000000
    -0.156539 0.938085 -0.309017
    -0.245485 0.969400 0.000000
    -0.452653 0.836429 -0.309017
    -0.546948 0.837166 0.000000
    -0.699715 0.644133 -0.309017
    -0.789141 0.614213 0.000000
    -0.870952 0.382035 -0.309017
    0.000000 -0.309017 -0.951056
    0.000000 0.000000 -1.000000
    0.100338 -0.292274 -0.951056
    0.189802 -0.243858 -0.951056
    0.258699 -0.169016 -0.951056
    0.299561 -0.075859 -0.951056
    0.307961 0.025518 -0.951056
    0.282989 0.124131 -0.951056
    0.227351 0.209292 -0.951056
    0.147076 0.271772 -0.951056
    0.050863 0.304802 -0.951056
    -0.050863 0.304802 -0.951056
    -0.147076 0.271772 -0.951056
    -0.227351 0.209292 -0.951056
    -0.282989 0.124131 -0.951056
    -0.307961 0.025518 -0.951056
    -0.299561 -0.075859 -0.951056
    -0.258699 -0.169016 -0.951056
    -0.189802 -0.243858 -0.951056
    -0.100338 -0.292274 -0.951056
    0.000000 -0.809017 -0.587785
    0.262687 -0.765182 -0.587785
    0.496909 -0.638428 -0.587785
    0.677282 -0.442490 -0.587785
    0.784261 -0.198602 -0.587785
    0.806254 0.066808 -0.587785
    0.740876 0.324978 -0.587785
    0.595213 0.547932 -0.587785
    0.385050 0.711509 -0.587785
    0.133160 0.797983 -0.587785
    -0.133160 0.797983 -0.587785
    -0.385050 0.711509 -0.587785
    -0.595213 0.547932 -0.587785
    -0.740876 0.324978 -0.587785
    -0.806254 0.066808 -0.587785
    -0.784261 -0.198602 -0.587785
    -0.677282 -0.442490 -0.587785
    -0.496909 -0.638428 -0.587785
    -0.262687 -0.765182 -0.587785
    0.000000 -0.587785 -0.809017
    0.190854 -0.555937 -0.809017
    0.361025 -0.463845 -0.809017
    0.492074 -0.321488 -0.809017
    0.569799 -0.144293 -0.809017
    0.585778 0.048539 -0.809017
    0.538278 0.236111 -0.809017
    0.432448 0.398096 -0.809017
    0.279755 0.516942 -0.809017
    0.096746 0.579769 -0.809017
    -0.096746 0.579769 -0.809017
    -0.279755 0.516942 -0.809017
    -0.432448 0.398096 -0.809017
    -0.538278 0.236111 -0.809017
    -0.585778 0.048539 -0.809017
    -0.569799 -0.144293 -0.809017
    -0.492074 -0.321488 -0.809017
    -0.361025 -0.463845 -0.809017
    -0.190854 -0.555937 -0.809017
    -1.000000 -0.000000 1.000000
    -0.945817 0.324700 1.000000
    -0.789141 0.614213 1.000000
    -0.546948 0.837166 1.000000
    -0.245485 0.969400 1.000000
    0.082579 0.996584 1.000000
    0.401695 0.915773 1.000000
    0.677282 0.735724 1.000000
    0.879474 0.475947 1.000000
    0.986361 0.164595 1.000000
    0.986361 -0.164595 1.000000
    0.879474 -0.475947 1.000000
    0.677282 -0.735724 1.000000
    0.401695 -0.915773 1.000000
    0.082579 -0.996584 1.000000
    -0.245485 -0.969400 1.000000
    -0.546948 -0.837166 1.000000
    -0.789141 -0.614213 1.000000
    -0.945817 -0.324700 1.000000
    -0.921954 -0.233471 1.309017
    -0.947808 0.078538 1.309017
    -0.870952 0.382035 1.309017
    -0.699715 0.644133 1.309017
    -0.452653 0.836429 1.309017
    -0.156539 0.938085 1.309017
    0.156539 0.938085 1.309017
    0.452653 0.836429 1.309017
    0.699715 0.644133 1.309017
    0.870952 0.382035 1.309017
    0.947808 0.078538 1.309017
    0.921954 -0.233471 1.309017
    0.796193 -0.520179 1.309017
    0.584151 -0.750517 1.309017
    0.308807 -0.899526 1.309017
    -0.000000 -0.951056 1.309017
    -0.308807 -0.899526 1.309017
    -0.584151 -0.750517 1.309017
    -0.796193 -0.520179 1.309017
    -0.000000 -0.309017 1.951057
    -0.000000 -0.000000 2.000000
    -0.100338 -0.292274 1.951057
    -0.189802 -0.243858 1.951057
    -0.258699 -0.169016 1.951057
    -0.299561 -0.075859 1.951057
    -0.307961 0.025518 1.951057
    -0.282989 0.124131 1.951057
    -0.227351 0.209292 1.951057
    -0.147076 0.271772 1.951057
    -0.050863 0.304802 1.951057
    0.050863 0.304802 1.951057
    0.147076 0.271772 1.951057
    0.227351 0.209292 1.951057
    0.282989 0.124131 1.951057
    0.307961 0.025518 1.951057
    0.299561 -0.075859 1.951057
    0.258699 -0.169016 1.951057
    0.189802 -0.243858 1.951057
    0.100338 -0.292274 1.951057
    -0.000000 -0.809017 1.587785
    -0.262687 -0.765182 1.587785
    -0.496909 -0.638428 1.587785
    -0.677282 -0.442490 1.587785
    -0.784261 -0.198602 1.587785
    -0.806254 0.066808 1.587785
    -0.740876 0.324978 1.587785
    -0.595213 0.547932 1.587785
    -0.385050 0.711509 1.587785
    -0.133160 0.797983 1.587785
    0.133160 0.797983 1.587785
    0.385050 0.711509 1.587785
    0.595213 0.547932 1.587785
    0.740876 0.324978 1.587785
    0.806254 0.066808 1.587785
    0.784261 -0.198602 1.587785
    0.677282 -0.442490 1.587785
    0.496909 -0.638428 1.587785
    0.262687 -0.765182 1.587785
    -0.000000 -0.587785 1.809017
    -0.190854 -0.555937 1.809017
    -0.361025 -0.463845 1.809017
    -0.492074 -0.321488 1.809017
    -0.569799 -0.144293 1.809017
    -0.585778 0.048539 1.809017
    -0.538278 0.236111 1.809017
    -0.432448 0.398096 1.809017
    -0.279755 0.516942 1.809017
    -0.096746 0.579769 1.809017
    0.096746 0.579769 1.809017
    0.279755 0.516942 1.809017
    0.432448 0.398096 1.809017
    0.538278 0.236111 1.809017
    0.585778 0.048539 1.809017
    0.569799 -0.144293 1.809017
    0.492074 -0.321488 1.809017
    0.361025 -0.463845 1.809017
    0.190854 -0.555937 1.809017
]


faces = [
    1 2 3
    3 2 4
    3 4 5
    5 4 6
    5 6 7
    7 6 8
    7 8 9
    9 8 10
    9 10 11
    11 10 12
    11 12 13
    13 12 14
    13 14 15
    15 14 16
    15 16 17
    17 16 18
    17 18 19
    19 18 20
    19 20 21
    21 20 22
    21 22 23
    23 22 24
    23 24 25
    25 24 26
    25 26 27
    27 26 28
    27 28 29
    29 28 30
    29 30 31
    31 30 32
    31 32 33
    33 32 34
    33 34 35
    35 34 36
    35 36 37
    37 36 38
    37 38 1
    1 38 2
    39 40 41
    41 40 42
    42 40 43
    43 40 44
    44 40 45
    45 40 46
    47 40 48
    48 40 49
    49 40 50
    50 40 51
    51 40 52
    53 40 54
    54 40 55
    55 40 56
    56 40 57
    57 40 58
    58 40 39
    12 59 14
    14 59 60
    14 60 16
    16 60 61
    16 61 18
    18 61 62
    18 62 20
    20 62 63
    20 63 22
    22 63 64
    22 64 24
    24 64 65
    24 65 26
    26 65 66
    26 66 28
    28 66 67
    28 67 30
    30 67 68
    30 68 32
    32 68 69
    32 69 34
    34 69 70
    34 70 36
    36 70 71
    36 71 38
    38 71 72
    38 72 2
    2 72 73
    2 73 4
    4 73 74
    4 74 6
    6 74 75
    6 75 8
    8 75 76
    8 76 10
    10 76 77
    10 77 12
    12 77 59
    59 78 60
    60 78 79
    60 79 61
    61 79 80
    61 80 62
    62 80 81
    62 81 63
    63 81 82
    63 82 64
    64 82 83
    64 83 65
    65 83 84
    65 84 66
    66 84 85
    66 85 67
    67 85 86
    67 86 68
    68 86 87
    68 87 69
    69 87 88
    69 88 70
    70 88 89
    70 89 71
    71 89 90
    71 90 72
    72 90 91
    72 91 73
    73 91 92
    73 92 74
    74 92 93
    74 93 75
    75 93 94
    75 94 76
    76 94 95
    76 95 77
    77 95 96
    77 96 59
    59 96 78
    78 39 79
    79 39 41
    79 41 80
    80 41 42
    80 42 81
    81 42 43
    81 43 82
    82 43 44
    82 44 83
    83 44 45
    83 45 84
    84 45 46
    84 46 85
    85 46 47
    85 47 86
    86 47 48
    86 48 87
    87 48 49
    87 49 88
    88 49 50
    88 50 89
    89 50 51
    89 51 90
    90 51 52
    90 52 91
    91 52 53
    91 53 92
    92 53 54
    92 54 93
    93 54 55
    93 55 94
    94 55 56
    94 56 95
    95 56 57
    95 57 96
    96 57 58
    96 58 78
    78 58 39
    46 40 47
    52 40 53
    5 97 3
    3 97 98
    3 98 1
    1 98 99
    1 99 37
    37 99 100
    37 100 35
    35 100 101
    35 101 33
    33 101 102
    33 102 31
    31 102 103
    31 103 29
    29 103 104
    29 104 27
    27 104 105
    27 105 25
    25 105 106
    25 106 23
    23 106 107
    23 107 21
    21 107 108
    21 108 19
    19 108 109
    19 109 17
    17 109 110
    17 110 15
    15 110 111
    15 111 13
    13 111 112
    13 112 11
    11 112 113
    11 113 9
    9 113 114
    9 114 7
    7 114 115
    7 115 5
    5 115 97
    115 116 97
    97 116 117
    97 117 98
    98 117 118
    98 118 99
    99 118 119
    99 119 100
    100 119 120
    100 120 101
    101 120 121
    101 121 102
    102 121 122
    102 122 103
    103 122 123
    103 123 104
    104 123 124
    104 124 105
    105 124 125
    105 125 106
    106 125 126
    106 126 107
    107 126 127
    107 127 108
    108 127 128
    108 128 109
    109 128 129
    109 129 110
    110 129 130
    110 130 111
    111 130 131
    111 131 112
    112 131 132
    112 132 113
    113 132 133
    113 133 114
    114 133 134
    114 134 115
    115 134 116
    135 136 137
    137 136 138
    138 136 139
    139 136 140
    140 136 141
    141 136 142
    142 136 143
    144 136 145
    145 136 146
    146 136 147
    147 136 148
    148 136 149
    150 136 151
    151 136 152
    152 136 153
    153 136 154
    154 136 135
    131 155 132
    132 155 156
    132 156 133
    133 156 157
    133 157 134
    134 157 158
    134 158 116
    116 158 159
    116 159 117
    117 159 160
    117 160 118
    118 160 161
    118 161 119
    119 161 162
    119 162 120
    120 162 163
    120 163 121
    121 163 164
    121 164 122
    122 164 165
    122 165 123
    123 165 166
    123 166 124
    124 166 167
    124 167 125
    125 167 168
    125 168 126
    126 168 169
    126 169 127
    127 169 170
    127 170 128
    128 170 171
    128 171 129
    129 171 172
    129 172 130
    130 172 173
    130 173 131
    131 173 155
    155 174 156
    156 174 175
    156 175 157
    157 175 176
    157 176 158
    158 176 177
    158 177 159
    159 177 178
    159 178 160
    160 178 179
    160 179 161
    161 179 180
    161 180 162
    162 180 181
    162 181 163
    163 181 182
    163 182 164
    164 182 183
    164 183 165
    165 183 184
    165 184 166
    166 184 185
    166 185 167
    167 185 186
    167 186 168
    168 186 187
    168 187 169
    169 187 188
    169 188 170
    170 188 189
    170 189 171
    171 189 190
    171 190 172
    172 190 191
    172 191 173
    173 191 192
    173 192 155
    155 192 174
    174 135 175
    175 135 137
    175 137 176
    176 137 138
    176 138 177
    177 138 139
    177 139 178
    178 139 140
    178 140 179
    179 140 141
    179 141 180
    180 141 142
    180 142 181
    181 142 143
    181 143 182
    182 143 144
    182 144 183
    183 144 145
    183 145 184
    184 145 146
    184 146 185
    185 146 147
    185 147 186
    186 147 148
    186 148 187
    187 148 149
    187 149 188
    188 149 150
    188 150 189
    189 150 151
    189 151 190
    190 151 152
    190 152 191
    191 152 153
    191 153 192
    192 153 154
    192 154 174
    174 154 135
    136 150 149
    136 144 143
]

# Index from 0, as expected by the C code
faces .-= 1


# Split vertices into the bottom and top caps
half = size(verts)[1] ÷ 2

verts_bottom = zeros(half, 3)
index_bottom = 1

verts_top = zeros(half, 3)
index_top = 1

faces_split = copy(faces)

for (i, v) in enumerate(eachrow(verts))
    if v[3] < 0.5
        verts_bottom[index_bottom, :] .= v
        faces_split[faces .== i - 1] .= index_bottom - 1

        global index_bottom += 1
    else
        verts_top[index_top, :] .= v
        faces_split[faces .== i - 1] .= index_top - 1 + half

        global index_top += 1
    end
end

verts_split = [verts_bottom ; verts_top]


# Move top cap back to (0, 0, 0) and flip it
verts_split[half + 1:end, 3] .-= 1
verts_split[half + 1:end, 3] .*= -1


# Write processed vertices and face indices to CSV files
open("spherical_cylinder_verts.csv", "w") do io
    writedlm(io, verts_split, ',')
end

open("spherical_cylinder_faces.csv", "w") do io
    writedlm(io, faces_split, ',')
end


# Create spherical cylinder spanning two points p1 and p2 on a particle trajectory
p1 = [5, 4, 3]
r1 = 2

p2 = [6, 7, 8]
r2 = 4


# Expand vertices to given particles' radii
verts_split[1:half, :] .*= r1
verts_split[half + 1:end, :] .*= r2


# Create rotation matrix from vertical +Z to |p2 - p1|
a = [0, 0, 1]
b = normalize(p2 .- p1)

v = cross(a, b)
c = dot(a, b)
vx = [  0    -v[3]   v[2]
       v[3]    0    -v[1]
      -v[2]   v[1]    0  ]

rot = I(3) + vx + vx * vx / (1 + c)

# Rotate all vertices
verts_split = (rot * verts_split')'


# Translate vertices
verts_split[1:half, :] .+= p1'
verts_split[half + 1:end, :] .+= p2'


# Plot mesh and direction from p1 to p2
m = GeometryBasics.Mesh(
    Point3.(eachrow(verts_split)),
    TriangleFace.(eachrow(faces_split .+ 1)),
)

wireframe(m)
lines!([p1 p2], color = :red)

current_figure()
