-
Notifications
You must be signed in to change notification settings - Fork 0
/
mueller.f90
265 lines (258 loc) · 9.07 KB
/
mueller.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
! Reads a csv from B23 Mueller polarimeter and performs Mueller matrix decomposition on each matrix in the array.
! There is a suspected bug in the matrix log algorithm due to mismatch with python program results.
! © Louis Minion 2023. Released under GNU GPL V3.
program mueller
implicit none
! Variable declaration
double precision :: G(4, 4) ! G = diag([1, -1, -1, -1]) Minkowski metric
double precision, allocatable, dimension( :, :, :) :: expmm, lu, lm, l
double precision, allocatable, dimension(:, :) :: expmmi, lui, lmi, li
external :: dgeev, DGETRF, DGETRI
intrinsic :: min, int, log, transpose, matmul
! double precision :: wavelengths, As, LDs,
Integer :: open_error, i,j, nlines, nargs
character(len=256) :: filename, outfilename
! Load matrix into memory
! Split matrix by wavelength
! For each wavelength, take the mueller matrix, find the logarithm
! Split the differential matrix into symmetric and asymmetric parts by mat multiplication
! Append an array for each of the variables extracted
type :: Experimental_Mueller_Matrix
integer :: record_n
double precision :: wavelength
double precision :: x_pos
double precision :: y_pos
character(19) :: time
double precision :: m_00
double precision :: m_01
double precision :: m_02
double precision :: m_03
double precision :: m_10
double precision :: m_11
double precision :: m_12
double precision :: m_13
double precision :: m_20
double precision :: m_21
double precision :: m_22
double precision :: m_23
double precision :: m_30
double precision :: m_31
double precision :: m_32
double precision :: m_33
double precision :: A
double precision :: LDabs
double precision :: LDang
double precision :: CD
double precision :: LRabs
double precision :: LRang
double precision :: CR
double precision :: LD
double precision :: LDp
double precision :: LR
double precision :: LRp
double precision :: TR
double precision :: AA
double precision :: LDabsA
double precision :: LDangA
double precision :: CDA
double precision :: LRabsA
double precision :: LRangA
double precision :: CRA
double precision :: LDA
double precision :: LDpA
double precision :: LRA
double precision :: LRpA
double precision :: TRA
double precision :: N
double precision :: NAngle
double precision :: S
double precision :: SAngle
double precision :: C
double precision :: Mraw00
double precision :: Mraw01
double precision :: Mraw02
double precision :: Mraw03
double precision :: Mraw10
double precision :: Mraw11
double precision :: Mraw12
double precision :: Mraw13
double precision :: Mraw20
double precision :: Mraw21
double precision :: Mraw22
double precision :: Mraw23
double precision :: Mraw30
double precision :: Mraw31
double precision :: Mraw32
double precision :: Mraw33
double precision :: MSTD00
double precision :: MSTD01
double precision :: MSTD02
double precision :: MSTD03
double precision :: MSTD10
double precision :: MSTD11
double precision :: MSTD12
double precision :: MSTD13
double precision :: MSTD20
double precision :: MSTD21
double precision :: MSTD22
double precision :: MSTD23
double precision :: MSTD30
double precision :: MSTD31
double precision :: MSTD32
double precision :: MSTD33
integer :: WavelengthIndex
integer :: XIndex
integer :: YIndex
integer :: TemperatureIndex
integer :: RepeatIndex
double precision :: Temp
end type Experimental_Mueller_Matrix
type(Experimental_Mueller_Matrix), allocatable, dimension(:) :: measured_matrices
nargs = command_argument_count()
if (nargs == 0) then
Write(*,*) 'No filename entered. Enter filename below: '
Read(*,*) filename
else
call GET_COMMAND_ARGUMENT(1,filename)
end if
filename = TRIM(filename)
j = LEN_TRIM(filename)
outfilename = filename(1:j-4) // '_LOGDECOMP.csv'
G = reshape((/double precision :: 1,0,0,0,0,-1,0,0,0,0,-1,0,0,0,0,-1/), (/4,4/))
nlines = 0
open(1, file = filename)
do
read(1,*, end=10)
nlines = nlines + 1
end do
10 close(1)
write(*,*) "Number of matrices in file:",nlines-17
allocate( measured_matrices(nlines-17))
Open(Unit=12, File=filename, Action='Read', Iostat=open_error, status='old')
If (open_error /= 0) Stop 'Error opening data file'
Write(*,*) 'Reading from input data...'
j = 0
do i= 1, 100000
if (j<17) then
read(12, *)
j = j +1
cycle
end if
j = j + 1
read(12, *, end=100) measured_matrices(i-17)
! write(*,*) measured_matrices(i-17)%wavelength
end do
100 Close(12)
allocate(expmm(4,4,nlines))
do i=1,nlines
expmm(1,1,i) = measured_matrices(i)%m_00
expmm(1,2,i) = measured_matrices(i)%m_01
expmm(1,3,i) = measured_matrices(i)%m_02
expmm(1,4,i) = measured_matrices(i)%m_03
expmm(2,1,i) = measured_matrices(i)%m_10
expmm(2,2,i) = measured_matrices(i)%m_11
expmm(2,3,i) = measured_matrices(i)%m_12
expmm(2,4,i) = measured_matrices(i)%m_13
expmm(3,1,i) = measured_matrices(i)%m_20
expmm(3,2,i) = measured_matrices(i)%m_21
expmm(3,3,i) = measured_matrices(i)%m_22
expmm(3,4,i) = measured_matrices(i)%m_23
expmm(4,1,i) = measured_matrices(i)%m_30
expmm(4,2,i) = measured_matrices(i)%m_31
expmm(4,3,i) = measured_matrices(i)%m_32
expmm(4,4,i) = measured_matrices(i)%m_33
end do
allocate(expmmi(4,4))
allocate(l(4,4,nlines-17))
allocate(li(4,4))
allocate(lui(4,4))
allocate(lmi(4,4))
allocate(lu(4,4,nlines-17))
allocate(lm(4,4,nlines-17))
do i=1, nlines-17
expmmi = expmm(:,:,i)
call CALC_LOGM(expmmi,li)
! write(*,*) li
l(:,:,i) = li
! Lm = 1/2(L-GL^TG), Lu = 1/2(L+GL^TG)
lmi = 0.5*(li-matmul(G,matmul(transpose(li),G)))
lui = 0.5*(li+matmul(G,matmul(transpose(li),G)))
lm(:,:,i) = lmi
lu(:,:,i) = lui
end do
! Write data out
open(unit=15,file=outfilename,Action='Write', Iostat=open_error)
If (open_error /= 0) Stop 'Error opening output file'
101 format(1x,*(g0, :, ", "))
do i=1,nlines-17
if (i==1) then
Write(15,101) "LM00","LM01","LM02","LM03","LM10","LM11", &
"LM12","LM13","LM20","LM21","LM22","LM23","LM30","LM31","LM32","LM33", &
"LU00","LU01","LU02","LU03","LU10","LU11","LU12","LU13","LU20","LU21", &
"LU22","LU23","LU30","LU31","LU32","LU33","L00","L01","L02","L03", &
"L10","L11","L12","L13","L20","L21","L22","L23","L30","L31","L32","L33"
end if
Write(15,101) lm(1,1,i) , lm(1,2,i), lm(1,3,i), lm(1,4,i), &
lm(2,1,i) , lm(2,2,i), lm(2,3,i), lm(2,4,i), &
lm(3,1,i) , lm(3,2,i), lm(3,3,i), lm(3,4,i), &
lm(4,1,i) , lm(4,2,i), lm(4,3,i), lm(4,4,i), &
lu(1,1,i) , lu(1,2,i), lu(1,3,i), lu(1,4,i), &
lu(2,1,i) , lu(2,2,i), lu(2,3,i), lu(2,4,i), &
lu(3,1,i) , lu(3,2,i), lu(3,3,i), lu(3,4,i), &
lu(4,1,i) , lu(4,2,i), lu(4,3,i), lu(4,4,i), &
l(1,1,i) , l(1,2,i), l(1,3,i), l(1,4,i), &
l(2,1,i) , l(2,2,i), l(2,3,i), l(2,4,i), &
l(3,1,i) , l(3,2,i), l(3,3,i), l(3,4,i), &
l(4,1,i) , l(4,2,i), l(4,3,i), l(4,4,i)
end do
Close(15)
end program mueller
SUBROUTINE CALC_LOGM(M, LOGM)
! M is an input square matrix to find the logarithm of. Must be diagonizable.
IMPLICIT NONE
integer :: n,i,j
Integer :: LWORK, LWMAX, INFO
Parameter (LWMAX = 1000)
Parameter (n=4)
double precision WORK(LWMAX)
double precision,dimension(n)::WR,WI,ipiv
double precision, dimension(n,n) :: M, LOGM, VL, VR, logADASH, VRINV
! external :: dgeev, dgetri, dgetrf
LWORK = -1
! OPTIMAL WORKSPACE QUERY
call dgeev('N','V',n,M,n, WR, WI, VL, n, VR, n,WORK,LWORK,INFO)
LWORK = MIN( LWMAX, INT( WORK( 1 ) ) )
call dgeev('N','V',n,M,n, WR, WI, VL, n, VR, n,WORK,LWORK,INFO)
if( INFO.GT.0 ) then
Write(*,*)'The algorithm failed to compute eigenvalues.'
GOTO 57
end if
! Check eigenvalues >0
do i=1,n
if (WR(i) .LE. 0.d0) then
Write(*,*) 'Eigenvalue less than zero, cannot compute '
GOTO 57
end if
end do
! initialise array as zero
do i=1,n
do j=1,n
logADASH(i,j) = 0.d0
end do
end do
do i =1,n
logADASH(i,i) = LOG(WR(i))
end do
VRINV = VR
call DGETRF(n,n,VRINV,n, ipiv, INFO)
if (info .NE. 0) then
Write(*,*) 'LU factorisation failed exit code', INFO
GOTO 57
end if
call DGETRI( n, VRINV, n, ipiv, WORK, LWORK, INFO )
if (info .NE. 0) then
Write(*,*) 'LU factorisation failed exit code', INFO
GOTO 57
end if
logM = matmul(matmul(VR,logADASH), VRINV)
57 END SUBROUTINE CALC_LOGM