← Back to branch summary

~ubuntu-branches/ubuntu/trusty/travis/trusty-proposed

« back to all changes in this revision

Viewing changes to src/normalmode.cpp

  • Committer: Package Import Robot
  • Author(s): Daniel Leidert
  • Date: 2014-01-18 20:07:16 UTC
  • mfrom: (1.1.8)
  • Revision ID: package-import@ubuntu.com-20140118200716-whsmcg7fa1eyqecq
Tags: 140117-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/normalmode.cpp
1
 
/*****************************************************************************
2
 
    TRAVIS - Trajectory Analyzer and Visualizer
3
 
    http://www.travis-analyzer.de/
4
 
 
5
 
    Copyright (c) 2009-2013 Martin Brehm
6
 
                  2012-2013 Martin Thomas
7
 
 
8
 
    This file written by Martin Thomas.
9
 
 
10
 
    This program is free software: you can redistribute it and/or modify
11
 
    it under the terms of the GNU General Public License as published by
12
 
    the Free Software Foundation, either version 3 of the License, or
13
 
    (at your option) any later version.
14
 
 
15
 
    This program is distributed in the hope that it will be useful,
16
 
    but WITHOUT ANY WARRANTY; without even the implied warranty of
17
 
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18
 
    GNU General Public License for more details.
19
 
 
20
 
    You should have received a copy of the GNU General Public License
21
 
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
22
 
*****************************************************************************/
23
 
 
24
 
#include "normalmode.h"
25
 
 
26
 
#include "acf.h"
27
 
#include "globalvar.h"
28
 
#include "tools.h"
29
 
#include "xfloatarray.h"
30
 
#include "xobarray.h"
31
 
 
32
 
#include <stdio.h>
33
 
 
34
 
static float g_freqStep;
35
 
static CxFloatArray *g_transMatrix;
36
 
 
37
 
static double findRotationAngle(int n, CxObArray *cc_matrix, int i, int j) {
38
 
        CxFloatArray *aii = (CxFloatArray *)cc_matrix->GetAt(i * n + i);
39
 
        CxFloatArray *aij = (CxFloatArray *)cc_matrix->GetAt(i * n + j);
40
 
        CxFloatArray *ajj = (CxFloatArray *)cc_matrix->GetAt(j * n + j);
41
 
        double int1 = 0.0, int2 = 0.0, int3 = 0.0;
42
 
        int k;
43
 
        for(k = 0; k < aij->GetSize(); k++) {
44
 
//              if(fabsf(aij->GetAt(k)) > 1.0e7f)
45
 
                        int1 += (double)aij->GetAt(k) * (double)aij->GetAt(k);
46
 
//              if((fabsf(aij->GetAt(k)) > 1.0e7f) && (fabsf(aii->GetAt(k)-ajj->GetAt(k)) > 1.0e7f))
47
 
                        int2 += (double)aij->GetAt(k) * ((double)aii->GetAt(k) - (double)ajj->GetAt(k));
48
 
//              if(fabsf(aii->GetAt(k)-ajj->GetAt(k)) > 1.0e7f)
49
 
                        int3 += ((double)aii->GetAt(k) - (double)ajj->GetAt(k)) * ((double)aii->GetAt(k) - (double)ajj->GetAt(k));
50
 
        }
51
 
        int1 *= (double)g_freqStep;
52
 
        int2 *= (double)g_freqStep;
53
 
        int3 *= (double)g_freqStep;
54
 
//      mprintf(RED, "---\n(%d, %d) int1 = %30g int2 = %30g int3 = %30g\n", i, j, int1, int2, int3);
55
 
        if(fabs(int2) < 1.0) {
56
 
//              mprintf(RED, "int2 < 1.0\n");
57
 
                double f = 4.0 * int1 - int3;
58
 
//              mprintf(RED, "f = %30g\n", f);
59
 
                if((fabs(int1) < 1.0) && (fabs(int2) < 1.0) && (fabs(int3) < 1.0))
60
 
                        return 0.0;
61
 
                if(0.5 * f > 0.0)
62
 
                        return 1.0;
63
 
                if(-2.0 * f > 0.0)
64
 
                        return 0.0;
65
 
        } else {
66
 
//              mprintf(RED, "int2 >= 1.0\n");
67
 
                double f = (4.0 * int1 - int3) / int2;
68
 
                double r[2][2];
69
 
                double s[2][2];
70
 
                for(k = 0; k < 2; k++) {
71
 
                        int l;
72
 
                        for(l = 0; l < 2; l++) {
73
 
//                              double r = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
74
 
//                              double s = -2.0 * int2 / pow(r * r + 1.0, 4.0) * (2.0 * r * (pow(r, 4.0) - 14.0 * r * r + 9.0) + f * (3.0 * pow(r, 4.0) - 8.0 * r * r + 1.0));
75
 
//                              if((fabs(r) <= 1.0) && (s > 0.0))
76
 
//                                      return (float)r;
77
 
                                r[k][l] = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
78
 
                                s[k][l] = -2.0 * int2 / pow(r[k][l] * r[k][l] + 1.0, 4.0) * (2.0 * r[k][l] * (pow(r[k][l], 4.0) - 14.0 * r[k][l] * r[k][l] + 9.0) + f * (3.0 * pow(r[k][l], 4.0) - 8.0 * r[k][l] * r[k][l] + 1.0));
79
 
//                              mprintf(RED, "t = %10g, s = %30g\n", r[k][l], s[k][l]);
80
 
                        }
81
 
                }
82
 
                for(k = 0; k < 2; k++) {
83
 
                        int l;
84
 
                        for(l = 0; l < 2; l++) {
85
 
                                if((fabs(r[k][l]) <= 1.0) && (s[k][l] > 0.0))
86
 
                                        return r[k][l];
87
 
                        }
88
 
                }
89
 
        }
90
 
//      double f = (4.0 * int1 - int3) / int2;
91
 
//      double r[2][2];
92
 
//      double s[2][2];
93
 
//      for(k = 0; k < 2; k++) {
94
 
//              int l;
95
 
//              for(l = 0; l < 2; l++) {
96
 
//                      double r = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
97
 
//                      double s = -2.0 * int2 / pow(r * r + 1.0, 4.0) * (2.0 * r * (pow(r, 4.0) - 14.0 * r * r + 9.0) + f * (3.0 * pow(r, 4.0) - 8.0 * r * r + 1.0));
98
 
//                      if((fabs(r) <= 1.0) && (s >= 0.0))
99
 
//                              return (float)r;
100
 
//                      r[k][l] = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
101
 
//                      s[k][l] = -2.0 * int2 / pow(r[k][l] * r[k][l] + 1.0, 4.0) * (2.0 * r[k][l] * (pow(r[k][l], 4.0) - 14.0 * r[k][l] * r[k][l] + 9.0) + f * (3.0 * pow(r[k][l], 4.0) - 8.0 * r[k][l] * r[k][l] + 1.0));
102
 
//                      mprintf(RED, "r = %10.6f, s = %30.2f\n", r[k][l], s[k][l]);
103
 
//              }
104
 
//      }
105
 
//      for(k = 0; k < 2; k++) {
106
 
//              int l;
107
 
//              for(l = 0; l < 2; l++) {
108
 
//                      if((fabs(r[k][l]) <= 1.0f) && (s[k][l] >= 0.0f))
109
 
//                              return r[k][l];
110
 
//              }
111
 
//      }
112
 
        
113
 
        mprintf(RED, "No suitable rotation angle found. Aborting.\n");
114
 
        return 0.1;
115
 
//      abort();
116
 
}
117
 
 
118
 
static double offDiagonalNorm(int n, CxObArray *cc_matrix) {
119
 
        double norm = 0.0;
120
 
        int i, j, k;
121
 
//      mprintf("@");
122
 
        for(i = 0; i < n; i++) {
123
 
                for(j = i + 1; j < n; j++) {
124
 
                        double integral = 0.0;
125
 
                        CxFloatArray *fa = (CxFloatArray *)cc_matrix->GetAt(i * n + j);
126
 
                        for(k = 0; k < fa->GetSize(); k++) {
127
 
//                              if(fabsf(fa->GetAt(k)) > 1.0e7f)
128
 
                                        integral += fa->GetAt(k) * fa->GetAt(k);
129
 
                        }
130
 
                        integral *= g_freqStep;
131
 
//                      mprintf("%14g ", integral);
132
 
                        norm += integral;
133
 
                }
134
 
//              mprintf("\n@");
135
 
        }
136
 
//      mprintf("\n");
137
 
        return sqrt(norm);
138
 
}
139
 
 
140
 
void normalModeAnalysis(int n, CxObArray *cc_matrix) {
141
 
        mprintf("    Fourier transforming cross correlation matrix...\n");
142
 
        mprintf(WHITE, "      [");
143
 
        int i, j, k;
144
 
 
145
 
        for(i = 0; i < n; i++) {
146
 
                for(j = 0; j < n; j++) {
147
 
                        CxFloatArray *temp;
148
 
                        try { temp = new CxFloatArray("normalModeAnalysis():temp"); } catch(...) { temp = NULL; }
149
 
                        if(temp == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
150
 
                        temp->CopyFrom((CxFloatArray *)cc_matrix->GetAt(i * n + j));
151
 
                        if(g_pGlobalVACF->m_bWindowFunction) {
152
 
                                for(k = 0; k < temp->GetSize(); k++) {
153
 
                                        temp->GetAt(k) *= powf(cosf((float)k / temp->GetSize() / 2.0f * Pi), 2.0f);
154
 
                                }
155
 
                        }
156
 
                        if(g_pGlobalVACF->m_iZeroPadding > 0) {
157
 
                                for(k = 0; k < g_pGlobalVACF->m_iZeroPadding; k++) {
158
 
                                        temp->Add(0.0f);
159
 
                                }
160
 
                        }
161
 
                        if(g_pGlobalVACF->m_iMirror != 0) {
162
 
                                int oldSize = temp->GetSize();
163
 
                                temp->SetSize(2 * oldSize);
164
 
                                for(k = 1; k < oldSize; k++) {
165
 
                                        temp->GetAt(oldSize + k) = temp->GetAt(oldSize - k);
166
 
                                }
167
 
                        }
168
 
                        CFFT *fft;
169
 
                        try { fft = new CFFT(); } catch(...) { fft = NULL; }
170
 
                        if(fft == NULL) NewException((double)sizeof(fft), __FILE__, __LINE__, __PRETTY_FUNCTION__);
171
 
                        fft->PrepareFFT_C2C(temp->GetSize());
172
 
                        for(k = 0; k < temp->GetSize(); k++) {
173
 
                                fft->m_pInput[2 * k] = temp->GetAt(k);
174
 
                                fft->m_pInput[2 * k + 1] = 0.0f;
175
 
                        }
176
 
                        fft->DoFFT();
177
 
                        ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->SetSize(temp->GetSize());
178
 
                        for(k = 0; k < temp->GetSize(); k++) {
179
 
                                ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = fft->m_pOutput[2 * k];
180
 
                        }
181
 
                        delete fft;
182
 
                        delete temp;
183
 
                        
184
 
//                      CxFloatArray *temp;
185
 
//                      try { temp = new CxFloatArray(); } catch(...) { temp = NULL; }
186
 
//                      if(temp == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
187
 
//                      temp->CopyFrom((CxFloatArray *)cc_matrix->GetAt(i * n + j));
188
 
//                      if(g_pGlobalVACF->m_bWindowFunction) {
189
 
//                              for(k = 0; k < temp->GetSize(); k++) {
190
 
//                                      temp->GetAt(k) *= powf(cosf((float)k / temp->GetSize() / 2.0f * Pi), 2.0f);
191
 
//                              }
192
 
//                      }
193
 
//                      if(g_pGlobalVACF->m_iZeroPadding > 0) {
194
 
//                              for(k = 0; k < g_pGlobalVACF->m_iZeroPadding; k++) {
195
 
//                                      temp->Add(0.0f);
196
 
//                              }
197
 
//                      }
198
 
//                      if(g_pGlobalVACF->m_iMirror != 0) {
199
 
//                              int oldSize = temp->GetSize();
200
 
//                              temp->SetSize(2 * oldSize);
201
 
//                              for(k = 1; k < oldSize; k++) {
202
 
//                                      temp->GetAt(oldSize + k) = temp->GetAt(oldSize - k);
203
 
//                              }
204
 
//                      }
205
 
//                      CFFT *fft;
206
 
//                      try { fft = new CFFT(); } catch(...) { fft = NULL; }
207
 
//                      if(fft == NULL) NewException((double)sizeof(fft), __FILE__, __LINE__, __PRETTY_FUNCTION__);
208
 
//                      fft->PrepareFFT_C2C(temp->GetSize());
209
 
//                      for(k = 0; k < temp->GetSize(); k++) {
210
 
//                              fft->m_pInput[2 * k] = temp->GetAt(k);
211
 
//                              fft->m_pInput[2 * k + 1] = 0.0f;
212
 
//                      }
213
 
//                      fft->DoFFT();
214
 
//                      ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->SetSize(temp->GetSize());
215
 
//                      for(k = 0; k < temp->GetSize(); k++) {
216
 
//                              ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = fft->m_pOutput[2 * k];
217
 
//                      }
218
 
//                      delete fft;
219
 
//                      delete temp;
220
 
//                      CACF *tempACF;
221
 
//                      try { tempACF = new CACF(); } catch(...) { tempACF = NULL; }
222
 
//                      if(tempACF == NULL) NewException((double)sizeof(CACF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
223
 
//                      tempACF->m_iSize = ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetSize();
224
 
//                      tempACF->m_bSpectrum = true;
225
 
//                      tempACF->m_bWindowFunction = g_pGlobalVACF->m_bWindowFunction;
226
 
//                      tempACF->m_fSpecWaveNumber = g_pGlobalVACF->m_fSpecWaveNumber;
227
 
//                      tempACF->m_iMirror = g_pGlobalVACF->m_iMirror;
228
 
//                      tempACF->m_iZeroPadding = g_pGlobalVACF->m_iZeroPadding;
229
 
//                      tempACF->Create();
230
 
//                      for(k = 0; k < tempACF->m_iSize; k++) {
231
 
//                              tempACF->m_pData[k] = ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k);
232
 
//                      }
233
 
//                      if(tempACF->m_iMirror != 0)
234
 
//                              tempACF->Mirror(tempACF->m_iMirror);
235
 
//                      if(tempACF->m_bWindowFunction)
236
 
//                              tempACF->Window();
237
 
//                      CFFT *fft;
238
 
//                      try { fft = new CFFT(); } catch(...) { fft = NULL; }
239
 
//                      if(fft == NULL) NewException((double)sizeof(fft), __FILE__, __LINE__, __PRETTY_FUNCTION__);
240
 
//                      fft->PrepareFFT_C2C(tempACF->m_iSize + tempACF->m_iZeroPadding);
241
 
//                      tempACF->Transform(fft);
242
 
//                      delete fft;
243
 
//                      ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->SetSize(tempACF->m_pSpectrum->m_iSize);
244
 
//                      for(k = 0; k < tempACF->m_pSpectrum->m_iSize; k++) {
245
 
//                              ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = tempACF->m_pSpectrum->m_pComplex[2*k];
246
 
//                      }
247
 
//                      delete tempACF;
248
 
                }
249
 
                mprintf(WHITE, "#");
250
 
        }
251
 
        mprintf(WHITE, "]\n");
252
 
        g_freqStep = 1e15f / 299792458.0f / 100.0f / g_fTimestepLength / g_iStride / ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize();
253
 
        
254
 
        mprintf("    Allocating transformation matrix...\n");
255
 
        try { g_transMatrix = new CxFloatArray("normalModeAnalysis():g_transMatrix"); } catch(...) { g_transMatrix = NULL; }
256
 
        if(g_transMatrix == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
257
 
//      g_transMatrix->SetSize(n * n);
258
 
        for(i = 0; i < n; i++) {
259
 
                for(j = 0; j < n; j++) {
260
 
                        g_transMatrix->Add((i == j) ? 1.0f : 0.0f);
261
 
                }
262
 
        }
263
 
        
264
 
        mprintf("    Saving transformation matrix\n");
265
 
        FILE *matrixi_file = OpenFileWrite("transformation_matrix_i.dat", false);
266
 
        for(i = 0; i < n; i++) {
267
 
                for(j = 0; j < n; j++) {
268
 
                        fprintf(matrixi_file, " %14g", g_transMatrix->GetAt(i * n + j));
269
 
                }
270
 
                fprintf(matrixi_file, "\n");
271
 
        }
272
 
        
273
 
        CxFloatArray *speciSum;
274
 
        try { speciSum = new CxFloatArray("normalModeAnalysis():speciSum"); } catch(...) { speciSum = NULL; }
275
 
        if(speciSum == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
276
 
        for(i = 0; i < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); i++)
277
 
                speciSum->Add(0.0f);
278
 
        for(i = 0; i < n; i++) {
279
 
                for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetSize(); k++) {
280
 
                        speciSum->GetAt(k) += ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k);
281
 
                }
282
 
        }
283
 
        FILE *speci_file = OpenFileWrite("normalmode_sum_initial.dat", false);
284
 
        for(k = 0;  k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
285
 
                fprintf(speci_file, "%14G %20G\n", k * g_freqStep, speciSum->GetAt(k));
286
 
        }
287
 
        fclose(speci_file);
288
 
        delete(speciSum);
289
 
        
290
 
        for(i = 0; i < n; i++) {
291
 
                for(j = i; j < n; j++) {
292
 
                        char name[128];
293
 
 
294
 
#ifdef TARGET_WINDOWS
295
 
                        _snprintf(name, 128, "speci_%d_%d.dat", i, j);
296
 
#elif defined TARGET_LINUX
297
 
                        snprintf(name, 128, "speci_%d_%d.dat", i, j);
298
 
#else
299
 
                        sprintf(name, "speci_%d_%d.dat", i, j);
300
 
#endif
301
 
 
302
 
                        FILE *speci_file = fopen(name, "w");
303
 
                        for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
304
 
                                fprintf(speci_file, "%14G %20G\n", k * g_freqStep, ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k));
305
 
                        }
306
 
                        fclose(speci_file);
307
 
                }
308
 
        }
309
 
        
310
 
        mprintf("    Minimizing cross correlations...\n");
311
 
        double norm = offDiagonalNorm(n, cc_matrix);
312
 
        mprintf("      %4d %20g\n", 0, norm);
313
 
        double change = norm;
314
 
        int count = 0;
315
 
        CxFloatArray *tempMatrix;
316
 
        try { tempMatrix = new CxFloatArray("normalModeAnalysis():tempMatrix"); } catch(...) { tempMatrix = NULL; }
317
 
        if(tempMatrix == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
318
 
        for(i = 0; i < n; i++) {
319
 
                tempMatrix->Add(0.0f);
320
 
        }
321
 
        while(fabs(change) > 1.0e-3 * norm) {
322
 
                count++;
323
 
                for(i = 0; i < n; i++) {
324
 
                        for(j = i + 1; j < n; j++) {
325
 
                                double t = findRotationAngle(n, cc_matrix, i, j);
326
 
                                double c = 1.0 / sqrt(t * t + 1.0);
327
 
                                double s = t * c;
328
 
//                              mprintf(RED, "(%d, %d) t = %10.6f c = %10.6f s = %10.6f\n", i, j, t, c, s);
329
 
                                for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); k++) {
330
 
                                        float temp[3];
331
 
                                        int l;
332
 
                                        for(l = 0; l < j; l++) {
333
 
                                                if(l < i) {
334
 
                                                        temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
335
 
                                                        temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
336
 
                                                        ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
337
 
                                                        ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
338
 
                                                } else if(l > i) {
339
 
                                                        tempMatrix->GetAt(l) = ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
340
 
                                                        temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
341
 
                                                        ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
342
 
                                                }
343
 
                                        }
344
 
                                        for(l = i + 1; l < n; l++) {
345
 
                                                if(l < j) {
346
 
                                                        temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * tempMatrix->GetAt(l);
347
 
                                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
348
 
                                                } else if(l > j) {
349
 
                                                        temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
350
 
                                                        temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
351
 
                                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
352
 
                                                        ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
353
 
                                                }
354
 
                                        }
355
 
                                        temp[0] = c * c * ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) + s * s * ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k) - 2 * c * s * ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k);
356
 
                                        temp[1] = (c * c - s * s) * ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) + c * s * (((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) - ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k));
357
 
                                        temp[2] = s * s * ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) + c * c * ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k) + 2 * c * s * ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k);
358
 
                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) = temp[0];
359
 
                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = temp[1];
360
 
                                        ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k) = temp[2];
361
 
                                        
362
 
//                                      for(l = 0; l < n; l++) {
363
 
//                                              if(l < i) {
364
 
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
365
 
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
366
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
367
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
368
 
//                                              } else if(l < j) {
369
 
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
370
 
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
371
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
372
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
373
 
//                                              } else {
374
 
////                                    for(l = 0; l < n; l++) {
375
 
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
376
 
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
377
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
378
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
379
 
//                                              }
380
 
//                                      }
381
 
//                                      for(l = n - 1; l >= 0; l--) {
382
 
//                                              if(l > j) {
383
 
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
384
 
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
385
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
386
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
387
 
//                                              } else if(l > i) {
388
 
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
389
 
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
390
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
391
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
392
 
//                                              } else {
393
 
////                                    for(l = 0; l < n; l++) {
394
 
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
395
 
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
396
 
////                                                    ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
397
 
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
398
 
//                                              }
399
 
//                                      }
400
 
                                }
401
 
                                float temp[2];
402
 
                                int l;
403
 
                                for(l = 0; l < n; l++) {
404
 
                                        temp[0] = c * g_transMatrix->GetAt(i * n + l) - s * g_transMatrix->GetAt(j * n + l);
405
 
                                        temp[1] = s * g_transMatrix->GetAt(i * n + l) + c * g_transMatrix->GetAt(j * n + l);
406
 
                                        g_transMatrix->GetAt(i * n + l) = temp[0];
407
 
                                        g_transMatrix->GetAt(j * n + l) = temp[1];
408
 
                                }
409
 
//                              mprintf("@@     (%d, %d) %20g\n", i, j, offDiagonalNorm(n, cc_matrix));
410
 
                        }
411
 
                }
412
 
                double newNorm = offDiagonalNorm(n, cc_matrix);
413
 
                change = newNorm - norm;
414
 
                norm = newNorm;
415
 
                mprintf("      %4d %20g\n", count, norm);
416
 
        }
417
 
        delete tempMatrix;
418
 
        
419
 
//      float temp[9];
420
 
//      for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); k++) {
421
 
//              temp[0] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
422
 
//              temp[1] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
423
 
//              temp[2] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
424
 
//              temp[3] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
425
 
//              temp[4] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
426
 
//              temp[5] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
427
 
//              temp[6] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
428
 
//              temp[7] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
429
 
//              temp[8] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
430
 
//              ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) = temp[0];
431
 
//              ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) = temp[1];
432
 
//              ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) = temp[2];
433
 
//              ((CxFloatArray *)cc_matrix->GetAt(12))->GetAt(k) = temp[4];
434
 
//              ((CxFloatArray *)cc_matrix->GetAt(15))->GetAt(k) = temp[5];
435
 
//              ((CxFloatArray *)cc_matrix->GetAt(24))->GetAt(k) = temp[8];
436
 
//      }
437
 
//      for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); k++) {
438
 
//              temp[0] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
439
 
//              temp[1] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
440
 
//              temp[2] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
441
 
//              temp[3] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
442
 
//              temp[4] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
443
 
//              temp[5] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
444
 
//              temp[6] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
445
 
//              temp[7] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
446
 
//              temp[8] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
447
 
//              ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) = temp[0];
448
 
//              ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) = temp[1];
449
 
//              ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) = temp[2];
450
 
//              ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) = temp[4];
451
 
//              ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) = temp[5];
452
 
//              ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k) = temp[8];
453
 
//      }
454
 
//      for(k = 0; k < 4096; k++) {
455
 
//              int l;
456
 
//              for (l=0;l<81;l++)
457
 
//              ((CxFloatArray *)cc_matrix->GetAt(l))->GetAt(k) = 1000.0f;
458
 
//      }
459
 
//      mprintf(RED, "Old: %20g, New: %20g\n", norm, offDiagonalNorm(n, cc_matrix));
460
 
        
461
 
        CxFloatArray *specSum;
462
 
        try { specSum = new CxFloatArray("normalModeAnalysis():specSum"); } catch(...) { specSum = NULL; }
463
 
        if(specSum == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
464
 
        for(i = 0; i < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); i++)
465
 
                specSum->Add(0.0f);
466
 
        
467
 
        mprintf("    Saving normal mode spectra...\n");
468
 
        for(i = 0; i < n; i++) {
469
 
                char name[128];
470
 
 
471
 
#ifdef TARGET_WINDOWS
472
 
                _snprintf(name, 128, "normalmode_%d.dat", i);
473
 
#elif defined TARGET_LINUX
474
 
                snprintf(name, 128, "normalmode_%d.dat", i);
475
 
#else
476
 
                sprintf(name, "normalmode_%d.dat", i);
477
 
#endif
478
 
 
479
 
                FILE *spec_file = OpenFileWrite(name, false);
480
 
                for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
481
 
                        fprintf(spec_file, "%14G %20G\n", k * g_freqStep, ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k));
482
 
                        specSum->GetAt(k) += ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k);
483
 
                }
484
 
                fclose(spec_file);
485
 
        }
486
 
        
487
 
        mprintf("    Saving sum of normal mode spectra...\n");
488
 
        FILE *spec_file = OpenFileWrite("normalmode_sum.dat", false);
489
 
        for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
490
 
                fprintf(spec_file, "%14G %20G\n", k * g_freqStep, specSum->GetAt(k));
491
 
        }
492
 
        fclose(spec_file);
493
 
        
494
 
        delete specSum;
495
 
        
496
 
        for(i = 0; i < n; i++) {
497
 
                for(j = i + 1; j < n; j++) {
498
 
                        char name[128];
499
 
 
500
 
#ifdef TARGET_WINDOWS
501
 
                        _snprintf(name, 128, "spec_%d_%d.dat", i, j);
502
 
#elif defined TARGET_LINUX
503
 
                        snprintf(name, 128, "spec_%d_%d.dat", i, j);
504
 
#else
505
 
                        sprintf(name, "spec_%d_%d.dat", i, j);
506
 
#endif
507
 
 
508
 
                        FILE *spec_file = fopen(name, "w");
509
 
                        for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
510
 
                                fprintf(spec_file, "%14G %20G\n", k * g_freqStep, ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k));
511
 
                        }
512
 
                        fclose(spec_file);
513
 
                }
514
 
        }
515
 
        
516
 
        mprintf("    Saving transformation matrix\n");
517
 
        FILE *matrix_file = OpenFileWrite("transformation_matrix.dat", false);
518
 
        for(i = 0; i < n; i++) {
519
 
                for(j = 0; j < n; j++) {
520
 
                        fprintf(matrix_file, " %14g", g_transMatrix->GetAt(i * n + j));
521
 
                }
522
 
                fprintf(matrix_file, "\n");
523
 
        }
524
 
        fclose(matrix_file);
525
 
        
526
 
        delete g_transMatrix;
527
 
}
 
1
/*****************************************************************************
 
2
    TRAVIS - Trajectory Analyzer and Visualizer
 
3
    http://www.travis-analyzer.de/
 
4
 
 
5
    Copyright (c) 2009-2014 Martin Brehm
 
6
                  2012-2014 Martin Thomas
 
7
 
 
8
    This file written by Martin Thomas.
 
9
 
 
10
    This program is free software: you can redistribute it and/or modify
 
11
    it under the terms of the GNU General Public License as published by
 
12
    the Free Software Foundation, either version 3 of the License, or
 
13
    (at your option) any later version.
 
14
 
 
15
    This program is distributed in the hope that it will be useful,
 
16
    but WITHOUT ANY WARRANTY; without even the implied warranty of
 
17
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
18
    GNU General Public License for more details.
 
19
 
 
20
    You should have received a copy of the GNU General Public License
 
21
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
22
*****************************************************************************/
 
23
 
 
24
#include "normalmode.h"
 
25
 
 
26
#include "acf.h"
 
27
#include "globalvar.h"
 
28
#include "tools.h"
 
29
#include "xfloatarray.h"
 
30
#include "xobarray.h"
 
31
 
 
32
#include <stdio.h>
 
33
 
 
34
static float g_freqStep;
 
35
static CxFloatArray *g_transMatrix;
 
36
 
 
37
static double findRotationAngle(int n, CxObArray *cc_matrix, int i, int j) {
 
38
        CxFloatArray *aii = (CxFloatArray *)cc_matrix->GetAt(i * n + i);
 
39
        CxFloatArray *aij = (CxFloatArray *)cc_matrix->GetAt(i * n + j);
 
40
        CxFloatArray *ajj = (CxFloatArray *)cc_matrix->GetAt(j * n + j);
 
41
        double int1 = 0.0, int2 = 0.0, int3 = 0.0;
 
42
        int k;
 
43
        for(k = 0; k < aij->GetSize(); k++) {
 
44
//              if(fabsf(aij->GetAt(k)) > 1.0e7f)
 
45
                        int1 += (double)aij->GetAt(k) * (double)aij->GetAt(k);
 
46
//              if((fabsf(aij->GetAt(k)) > 1.0e7f) && (fabsf(aii->GetAt(k)-ajj->GetAt(k)) > 1.0e7f))
 
47
                        int2 += (double)aij->GetAt(k) * ((double)aii->GetAt(k) - (double)ajj->GetAt(k));
 
48
//              if(fabsf(aii->GetAt(k)-ajj->GetAt(k)) > 1.0e7f)
 
49
                        int3 += ((double)aii->GetAt(k) - (double)ajj->GetAt(k)) * ((double)aii->GetAt(k) - (double)ajj->GetAt(k));
 
50
        }
 
51
        int1 *= (double)g_freqStep;
 
52
        int2 *= (double)g_freqStep;
 
53
        int3 *= (double)g_freqStep;
 
54
//      mprintf(RED, "---\n(%d, %d) int1 = %30g int2 = %30g int3 = %30g\n", i, j, int1, int2, int3);
 
55
        if(fabs(int2) < 1.0) {
 
56
//              mprintf(RED, "int2 < 1.0\n");
 
57
                double f = 4.0 * int1 - int3;
 
58
//              mprintf(RED, "f = %30g\n", f);
 
59
                if((fabs(int1) < 1.0) && (fabs(int2) < 1.0) && (fabs(int3) < 1.0))
 
60
                        return 0.0;
 
61
                if(0.5 * f > 0.0)
 
62
                        return 1.0;
 
63
                if(-2.0 * f > 0.0)
 
64
                        return 0.0;
 
65
        } else {
 
66
//              mprintf(RED, "int2 >= 1.0\n");
 
67
                double f = (4.0 * int1 - int3) / int2;
 
68
                double r[2][2];
 
69
                double s[2][2];
 
70
                for(k = 0; k < 2; k++) {
 
71
                        int l;
 
72
                        for(l = 0; l < 2; l++) {
 
73
//                              double r = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
 
74
//                              double s = -2.0 * int2 / pow(r * r + 1.0, 4.0) * (2.0 * r * (pow(r, 4.0) - 14.0 * r * r + 9.0) + f * (3.0 * pow(r, 4.0) - 8.0 * r * r + 1.0));
 
75
//                              if((fabs(r) <= 1.0) && (s > 0.0))
 
76
//                                      return (float)r;
 
77
                                r[k][l] = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
 
78
                                s[k][l] = -2.0 * int2 / pow(r[k][l] * r[k][l] + 1.0, 4.0) * (2.0 * r[k][l] * (pow(r[k][l], 4.0) - 14.0 * r[k][l] * r[k][l] + 9.0) + f * (3.0 * pow(r[k][l], 4.0) - 8.0 * r[k][l] * r[k][l] + 1.0));
 
79
//                              mprintf(RED, "t = %10g, s = %30g\n", r[k][l], s[k][l]);
 
80
                        }
 
81
                }
 
82
                for(k = 0; k < 2; k++) {
 
83
                        int l;
 
84
                        for(l = 0; l < 2; l++) {
 
85
                                if((fabs(r[k][l]) <= 1.0) && (s[k][l] > 0.0))
 
86
                                        return r[k][l];
 
87
                        }
 
88
                }
 
89
        }
 
90
//      double f = (4.0 * int1 - int3) / int2;
 
91
//      double r[2][2];
 
92
//      double s[2][2];
 
93
//      for(k = 0; k < 2; k++) {
 
94
//              int l;
 
95
//              for(l = 0; l < 2; l++) {
 
96
//                      double r = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
 
97
//                      double s = -2.0 * int2 / pow(r * r + 1.0, 4.0) * (2.0 * r * (pow(r, 4.0) - 14.0 * r * r + 9.0) + f * (3.0 * pow(r, 4.0) - 8.0 * r * r + 1.0));
 
98
//                      if((fabs(r) <= 1.0) && (s >= 0.0))
 
99
//                              return (float)r;
 
100
//                      r[k][l] = (-f + (k == 0 ? -1.0 : 1.0) * sqrt(16.0 + f * f) + (l == 0 ? -1.0 : 1.0) * sqrt(2.0) * sqrt(16.0 + f * f - (k == 0 ? -1.0 : 1.0) * f * sqrt(16.0 + f * f))) / 4.0;
 
101
//                      s[k][l] = -2.0 * int2 / pow(r[k][l] * r[k][l] + 1.0, 4.0) * (2.0 * r[k][l] * (pow(r[k][l], 4.0) - 14.0 * r[k][l] * r[k][l] + 9.0) + f * (3.0 * pow(r[k][l], 4.0) - 8.0 * r[k][l] * r[k][l] + 1.0));
 
102
//                      mprintf(RED, "r = %10.6f, s = %30.2f\n", r[k][l], s[k][l]);
 
103
//              }
 
104
//      }
 
105
//      for(k = 0; k < 2; k++) {
 
106
//              int l;
 
107
//              for(l = 0; l < 2; l++) {
 
108
//                      if((fabs(r[k][l]) <= 1.0f) && (s[k][l] >= 0.0f))
 
109
//                              return r[k][l];
 
110
//              }
 
111
//      }
 
112
        
 
113
        mprintf(RED, "No suitable rotation angle found. Aborting.\n");
 
114
        return 0.1;
 
115
//      abort();
 
116
}
 
117
 
 
118
static double offDiagonalNorm(int n, CxObArray *cc_matrix) {
 
119
        double norm = 0.0;
 
120
        int i, j, k;
 
121
//      mprintf("@");
 
122
        for(i = 0; i < n; i++) {
 
123
                for(j = i + 1; j < n; j++) {
 
124
                        double integral = 0.0;
 
125
                        CxFloatArray *fa = (CxFloatArray *)cc_matrix->GetAt(i * n + j);
 
126
                        for(k = 0; k < fa->GetSize(); k++) {
 
127
//                              if(fabsf(fa->GetAt(k)) > 1.0e7f)
 
128
                                        integral += fa->GetAt(k) * fa->GetAt(k);
 
129
                        }
 
130
                        integral *= g_freqStep;
 
131
//                      mprintf("%14g ", integral);
 
132
                        norm += integral;
 
133
                }
 
134
//              mprintf("\n@");
 
135
        }
 
136
//      mprintf("\n");
 
137
        return sqrt(norm);
 
138
}
 
139
 
 
140
void normalModeAnalysis(int n, CxObArray *cc_matrix) {
 
141
        mprintf("    Fourier transforming cross correlation matrix...\n");
 
142
        mprintf(WHITE, "      [");
 
143
        int i, j, k;
 
144
 
 
145
        for(i = 0; i < n; i++) {
 
146
                for(j = 0; j < n; j++) {
 
147
                        CxFloatArray *temp;
 
148
                        try { temp = new CxFloatArray("normalModeAnalysis():temp"); } catch(...) { temp = NULL; }
 
149
                        if(temp == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
150
                        temp->CopyFrom((CxFloatArray *)cc_matrix->GetAt(i * n + j));
 
151
                        if(g_pGlobalVACF->m_bWindowFunction) {
 
152
                                for(k = 0; k < temp->GetSize(); k++) {
 
153
                                        temp->GetAt(k) *= powf(cosf((float)k / temp->GetSize() / 2.0f * Pi), 2.0f);
 
154
                                }
 
155
                        }
 
156
                        if(g_pGlobalVACF->m_iZeroPadding > 0) {
 
157
                                for(k = 0; k < g_pGlobalVACF->m_iZeroPadding; k++) {
 
158
                                        temp->Add(0.0f);
 
159
                                }
 
160
                        }
 
161
                        if(g_pGlobalVACF->m_iMirror != 0) {
 
162
                                int oldSize = temp->GetSize();
 
163
                                temp->SetSize(2 * oldSize);
 
164
                                for(k = 1; k < oldSize; k++) {
 
165
                                        temp->GetAt(oldSize + k) = temp->GetAt(oldSize - k);
 
166
                                }
 
167
                        }
 
168
                        CFFT *fft;
 
169
                        try { fft = new CFFT(); } catch(...) { fft = NULL; }
 
170
                        if(fft == NULL) NewException((double)sizeof(fft), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
171
                        fft->PrepareFFT_C2C(temp->GetSize());
 
172
                        for(k = 0; k < temp->GetSize(); k++) {
 
173
                                fft->m_pInput[2 * k] = temp->GetAt(k);
 
174
                                fft->m_pInput[2 * k + 1] = 0.0f;
 
175
                        }
 
176
                        fft->DoFFT();
 
177
                        ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->SetSize(temp->GetSize());
 
178
                        for(k = 0; k < temp->GetSize(); k++) {
 
179
                                ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = fft->m_pOutput[2 * k];
 
180
                        }
 
181
                        delete fft;
 
182
                        delete temp;
 
183
                        
 
184
//                      CxFloatArray *temp;
 
185
//                      try { temp = new CxFloatArray(); } catch(...) { temp = NULL; }
 
186
//                      if(temp == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
187
//                      temp->CopyFrom((CxFloatArray *)cc_matrix->GetAt(i * n + j));
 
188
//                      if(g_pGlobalVACF->m_bWindowFunction) {
 
189
//                              for(k = 0; k < temp->GetSize(); k++) {
 
190
//                                      temp->GetAt(k) *= powf(cosf((float)k / temp->GetSize() / 2.0f * Pi), 2.0f);
 
191
//                              }
 
192
//                      }
 
193
//                      if(g_pGlobalVACF->m_iZeroPadding > 0) {
 
194
//                              for(k = 0; k < g_pGlobalVACF->m_iZeroPadding; k++) {
 
195
//                                      temp->Add(0.0f);
 
196
//                              }
 
197
//                      }
 
198
//                      if(g_pGlobalVACF->m_iMirror != 0) {
 
199
//                              int oldSize = temp->GetSize();
 
200
//                              temp->SetSize(2 * oldSize);
 
201
//                              for(k = 1; k < oldSize; k++) {
 
202
//                                      temp->GetAt(oldSize + k) = temp->GetAt(oldSize - k);
 
203
//                              }
 
204
//                      }
 
205
//                      CFFT *fft;
 
206
//                      try { fft = new CFFT(); } catch(...) { fft = NULL; }
 
207
//                      if(fft == NULL) NewException((double)sizeof(fft), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
208
//                      fft->PrepareFFT_C2C(temp->GetSize());
 
209
//                      for(k = 0; k < temp->GetSize(); k++) {
 
210
//                              fft->m_pInput[2 * k] = temp->GetAt(k);
 
211
//                              fft->m_pInput[2 * k + 1] = 0.0f;
 
212
//                      }
 
213
//                      fft->DoFFT();
 
214
//                      ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->SetSize(temp->GetSize());
 
215
//                      for(k = 0; k < temp->GetSize(); k++) {
 
216
//                              ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = fft->m_pOutput[2 * k];
 
217
//                      }
 
218
//                      delete fft;
 
219
//                      delete temp;
 
220
//                      CACF *tempACF;
 
221
//                      try { tempACF = new CACF(); } catch(...) { tempACF = NULL; }
 
222
//                      if(tempACF == NULL) NewException((double)sizeof(CACF), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
223
//                      tempACF->m_iSize = ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetSize();
 
224
//                      tempACF->m_bSpectrum = true;
 
225
//                      tempACF->m_bWindowFunction = g_pGlobalVACF->m_bWindowFunction;
 
226
//                      tempACF->m_fSpecWaveNumber = g_pGlobalVACF->m_fSpecWaveNumber;
 
227
//                      tempACF->m_iMirror = g_pGlobalVACF->m_iMirror;
 
228
//                      tempACF->m_iZeroPadding = g_pGlobalVACF->m_iZeroPadding;
 
229
//                      tempACF->Create();
 
230
//                      for(k = 0; k < tempACF->m_iSize; k++) {
 
231
//                              tempACF->m_pData[k] = ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k);
 
232
//                      }
 
233
//                      if(tempACF->m_iMirror != 0)
 
234
//                              tempACF->Mirror(tempACF->m_iMirror);
 
235
//                      if(tempACF->m_bWindowFunction)
 
236
//                              tempACF->Window();
 
237
//                      CFFT *fft;
 
238
//                      try { fft = new CFFT(); } catch(...) { fft = NULL; }
 
239
//                      if(fft == NULL) NewException((double)sizeof(fft), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
240
//                      fft->PrepareFFT_C2C(tempACF->m_iSize + tempACF->m_iZeroPadding);
 
241
//                      tempACF->Transform(fft);
 
242
//                      delete fft;
 
243
//                      ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->SetSize(tempACF->m_pSpectrum->m_iSize);
 
244
//                      for(k = 0; k < tempACF->m_pSpectrum->m_iSize; k++) {
 
245
//                              ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = tempACF->m_pSpectrum->m_pComplex[2*k];
 
246
//                      }
 
247
//                      delete tempACF;
 
248
                }
 
249
                mprintf(WHITE, "#");
 
250
        }
 
251
        mprintf(WHITE, "]\n");
 
252
        g_freqStep = 1e15f / 299792458.0f / 100.0f / g_fTimestepLength / g_iStride / ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize();
 
253
        
 
254
        mprintf("    Allocating transformation matrix...\n");
 
255
        try { g_transMatrix = new CxFloatArray("normalModeAnalysis():g_transMatrix"); } catch(...) { g_transMatrix = NULL; }
 
256
        if(g_transMatrix == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
257
//      g_transMatrix->SetSize(n * n);
 
258
        for(i = 0; i < n; i++) {
 
259
                for(j = 0; j < n; j++) {
 
260
                        g_transMatrix->Add((i == j) ? 1.0f : 0.0f);
 
261
                }
 
262
        }
 
263
        
 
264
        mprintf("    Saving transformation matrix\n");
 
265
        FILE *matrixi_file = OpenFileWrite("transformation_matrix_i.dat", false);
 
266
        for(i = 0; i < n; i++) {
 
267
                for(j = 0; j < n; j++) {
 
268
                        fprintf(matrixi_file, " %14g", g_transMatrix->GetAt(i * n + j));
 
269
                }
 
270
                fprintf(matrixi_file, "\n");
 
271
        }
 
272
        
 
273
        CxFloatArray *speciSum;
 
274
        try { speciSum = new CxFloatArray("normalModeAnalysis():speciSum"); } catch(...) { speciSum = NULL; }
 
275
        if(speciSum == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
276
        for(i = 0; i < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); i++)
 
277
                speciSum->Add(0.0f);
 
278
        for(i = 0; i < n; i++) {
 
279
                for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetSize(); k++) {
 
280
                        speciSum->GetAt(k) += ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k);
 
281
                }
 
282
        }
 
283
        FILE *speci_file = OpenFileWrite("normalmode_sum_initial.dat", false);
 
284
        for(k = 0;  k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
 
285
                fprintf(speci_file, "%14G %20G\n", k * g_freqStep, speciSum->GetAt(k));
 
286
        }
 
287
        fclose(speci_file);
 
288
        delete(speciSum);
 
289
        
 
290
        for(i = 0; i < n; i++) {
 
291
                for(j = i; j < n; j++) {
 
292
                        char name[128];
 
293
 
 
294
#ifdef TARGET_WINDOWS
 
295
                        _snprintf(name, 128, "speci_%d_%d.dat", i, j);
 
296
#elif defined TARGET_LINUX
 
297
                        snprintf(name, 128, "speci_%d_%d.dat", i, j);
 
298
#else
 
299
                        sprintf(name, "speci_%d_%d.dat", i, j);
 
300
#endif
 
301
 
 
302
                        FILE *speci_file = fopen(name, "w");
 
303
                        for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
 
304
                                fprintf(speci_file, "%14G %20G\n", k * g_freqStep, ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k));
 
305
                        }
 
306
                        fclose(speci_file);
 
307
                }
 
308
        }
 
309
        
 
310
        mprintf("    Minimizing cross correlations...\n");
 
311
        double norm = offDiagonalNorm(n, cc_matrix);
 
312
        mprintf("      %4d %20g\n", 0, norm);
 
313
        double change = norm;
 
314
        int count = 0;
 
315
        CxFloatArray *tempMatrix;
 
316
        try { tempMatrix = new CxFloatArray("normalModeAnalysis():tempMatrix"); } catch(...) { tempMatrix = NULL; }
 
317
        if(tempMatrix == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
318
        for(i = 0; i < n; i++) {
 
319
                tempMatrix->Add(0.0f);
 
320
        }
 
321
        while(fabs(change) > 1.0e-3 * norm) {
 
322
                count++;
 
323
                for(i = 0; i < n; i++) {
 
324
                        for(j = i + 1; j < n; j++) {
 
325
                                double t = findRotationAngle(n, cc_matrix, i, j);
 
326
                                double c = 1.0 / sqrt(t * t + 1.0);
 
327
                                double s = t * c;
 
328
//                              mprintf(RED, "(%d, %d) t = %10.6f c = %10.6f s = %10.6f\n", i, j, t, c, s);
 
329
                                for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); k++) {
 
330
                                        float temp[3];
 
331
                                        int l;
 
332
                                        for(l = 0; l < j; l++) {
 
333
                                                if(l < i) {
 
334
                                                        temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
335
                                                        temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
336
                                                        ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
 
337
                                                        ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
 
338
                                                } else if(l > i) {
 
339
                                                        tempMatrix->GetAt(l) = ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
340
                                                        temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
341
                                                        ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
 
342
                                                }
 
343
                                        }
 
344
                                        for(l = i + 1; l < n; l++) {
 
345
                                                if(l < j) {
 
346
                                                        temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * tempMatrix->GetAt(l);
 
347
                                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
 
348
                                                } else if(l > j) {
 
349
                                                        temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
 
350
                                                        temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
 
351
                                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
 
352
                                                        ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
 
353
                                                }
 
354
                                        }
 
355
                                        temp[0] = c * c * ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) + s * s * ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k) - 2 * c * s * ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k);
 
356
                                        temp[1] = (c * c - s * s) * ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) + c * s * (((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) - ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k));
 
357
                                        temp[2] = s * s * ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) + c * c * ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k) + 2 * c * s * ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k);
 
358
                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k) = temp[0];
 
359
                                        ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k) = temp[1];
 
360
                                        ((CxFloatArray *)cc_matrix->GetAt(j * n + j))->GetAt(k) = temp[2];
 
361
                                        
 
362
//                                      for(l = 0; l < n; l++) {
 
363
//                                              if(l < i) {
 
364
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
365
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
366
//                                                      ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
 
367
//                                                      ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
 
368
//                                              } else if(l < j) {
 
369
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
370
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
371
//                                                      ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
 
372
//                                                      ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
 
373
//                                              } else {
 
374
////                                    for(l = 0; l < n; l++) {
 
375
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
 
376
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
 
377
//                                                      ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) = temp[0];
 
378
//                                                      ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k) = temp[1];
 
379
//                                              }
 
380
//                                      }
 
381
//                                      for(l = n - 1; l >= 0; l--) {
 
382
//                                              if(l > j) {
 
383
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
 
384
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(j * n + l))->GetAt(k);
 
385
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
 
386
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
 
387
//                                              } else if(l > i) {
 
388
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
389
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(i * n + l))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
390
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
 
391
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
 
392
//                                              } else {
 
393
////                                    for(l = 0; l < n; l++) {
 
394
//                                                      temp[0] = c * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) - s * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
395
//                                                      temp[1] = s * ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) + c * ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k);
 
396
////                                                    ((CxFloatArray *)cc_matrix->GetAt(l * n + i))->GetAt(k) = temp[0];
 
397
//                                                      ((CxFloatArray *)cc_matrix->GetAt(l * n + j))->GetAt(k) = temp[1];
 
398
//                                              }
 
399
//                                      }
 
400
                                }
 
401
                                float temp[2];
 
402
                                int l;
 
403
                                for(l = 0; l < n; l++) {
 
404
                                        temp[0] = c * g_transMatrix->GetAt(i * n + l) - s * g_transMatrix->GetAt(j * n + l);
 
405
                                        temp[1] = s * g_transMatrix->GetAt(i * n + l) + c * g_transMatrix->GetAt(j * n + l);
 
406
                                        g_transMatrix->GetAt(i * n + l) = temp[0];
 
407
                                        g_transMatrix->GetAt(j * n + l) = temp[1];
 
408
                                }
 
409
//                              mprintf("@@     (%d, %d) %20g\n", i, j, offDiagonalNorm(n, cc_matrix));
 
410
                        }
 
411
                }
 
412
                double newNorm = offDiagonalNorm(n, cc_matrix);
 
413
                change = newNorm - norm;
 
414
                norm = newNorm;
 
415
                mprintf("      %4d %20g\n", count, norm);
 
416
        }
 
417
        delete tempMatrix;
 
418
        
 
419
//      float temp[9];
 
420
//      for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); k++) {
 
421
//              temp[0] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
 
422
//              temp[1] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
 
423
//              temp[2] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
 
424
//              temp[3] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
 
425
//              temp[4] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
 
426
//              temp[5] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
 
427
//              temp[6] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
 
428
//              temp[7] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
 
429
//              temp[8] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
 
430
//              ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) = temp[0];
 
431
//              ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) = temp[1];
 
432
//              ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) = temp[2];
 
433
//              ((CxFloatArray *)cc_matrix->GetAt(12))->GetAt(k) = temp[4];
 
434
//              ((CxFloatArray *)cc_matrix->GetAt(15))->GetAt(k) = temp[5];
 
435
//              ((CxFloatArray *)cc_matrix->GetAt(24))->GetAt(k) = temp[8];
 
436
//      }
 
437
//      for(k = 0; k < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); k++) {
 
438
//              temp[0] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
 
439
//              temp[1] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
 
440
//              temp[2] = -0.0259f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) - 0.978f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.209f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
 
441
//              temp[3] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
 
442
//              temp[4] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
 
443
//              temp[5] = 0.986f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.009f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.164f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
 
444
//              temp[6] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k);
 
445
//              temp[7] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k);
 
446
//              temp[8] = -0.163f * ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) + 0.210f * ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) + 0.964f * ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k);
 
447
//              ((CxFloatArray *)cc_matrix->GetAt(0))->GetAt(k) = temp[0];
 
448
//              ((CxFloatArray *)cc_matrix->GetAt(3))->GetAt(k) = temp[1];
 
449
//              ((CxFloatArray *)cc_matrix->GetAt(6))->GetAt(k) = temp[2];
 
450
//              ((CxFloatArray *)cc_matrix->GetAt(30))->GetAt(k) = temp[4];
 
451
//              ((CxFloatArray *)cc_matrix->GetAt(33))->GetAt(k) = temp[5];
 
452
//              ((CxFloatArray *)cc_matrix->GetAt(60))->GetAt(k) = temp[8];
 
453
//      }
 
454
//      for(k = 0; k < 4096; k++) {
 
455
//              int l;
 
456
//              for (l=0;l<81;l++)
 
457
//              ((CxFloatArray *)cc_matrix->GetAt(l))->GetAt(k) = 1000.0f;
 
458
//      }
 
459
//      mprintf(RED, "Old: %20g, New: %20g\n", norm, offDiagonalNorm(n, cc_matrix));
 
460
        
 
461
        CxFloatArray *specSum;
 
462
        try { specSum = new CxFloatArray("normalModeAnalysis():specSum"); } catch(...) { specSum = NULL; }
 
463
        if(specSum == NULL) NewException((double)sizeof(CxFloatArray), __FILE__, __LINE__, __PRETTY_FUNCTION__);
 
464
        for(i = 0; i < ((CxFloatArray *)cc_matrix->GetAt(0))->GetSize(); i++)
 
465
                specSum->Add(0.0f);
 
466
        
 
467
        mprintf("    Saving normal mode spectra...\n");
 
468
        for(i = 0; i < n; i++) {
 
469
                char name[128];
 
470
 
 
471
#ifdef TARGET_WINDOWS
 
472
                _snprintf(name, 128, "normalmode_%d.dat", i);
 
473
#elif defined TARGET_LINUX
 
474
                snprintf(name, 128, "normalmode_%d.dat", i);
 
475
#else
 
476
                sprintf(name, "normalmode_%d.dat", i);
 
477
#endif
 
478
 
 
479
                FILE *spec_file = OpenFileWrite(name, false);
 
480
                for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
 
481
                        fprintf(spec_file, "%14G %20G\n", k * g_freqStep, ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k));
 
482
                        specSum->GetAt(k) += ((CxFloatArray *)cc_matrix->GetAt(i * n + i))->GetAt(k);
 
483
                }
 
484
                fclose(spec_file);
 
485
        }
 
486
        
 
487
        mprintf("    Saving sum of normal mode spectra...\n");
 
488
        FILE *spec_file = OpenFileWrite("normalmode_sum.dat", false);
 
489
        for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
 
490
                fprintf(spec_file, "%14G %20G\n", k * g_freqStep, specSum->GetAt(k));
 
491
        }
 
492
        fclose(spec_file);
 
493
        
 
494
        delete specSum;
 
495
        
 
496
        for(i = 0; i < n; i++) {
 
497
                for(j = i + 1; j < n; j++) {
 
498
                        char name[128];
 
499
 
 
500
#ifdef TARGET_WINDOWS
 
501
                        _snprintf(name, 128, "spec_%d_%d.dat", i, j);
 
502
#elif defined TARGET_LINUX
 
503
                        snprintf(name, 128, "spec_%d_%d.dat", i, j);
 
504
#else
 
505
                        sprintf(name, "spec_%d_%d.dat", i, j);
 
506
#endif
 
507
 
 
508
                        FILE *spec_file = fopen(name, "w");
 
509
                        for(k = 0; k * g_freqStep < g_pGlobalVACF->m_fSpecWaveNumber; k++) {
 
510
                                fprintf(spec_file, "%14G %20G\n", k * g_freqStep, ((CxFloatArray *)cc_matrix->GetAt(i * n + j))->GetAt(k));
 
511
                        }
 
512
                        fclose(spec_file);
 
513
                }
 
514
        }
 
515
        
 
516
        mprintf("    Saving transformation matrix\n");
 
517
        FILE *matrix_file = OpenFileWrite("transformation_matrix.dat", false);
 
518
        for(i = 0; i < n; i++) {
 
519
                for(j = 0; j < n; j++) {
 
520
                        fprintf(matrix_file, " %14g", g_transMatrix->GetAt(i * n + j));
 
521
                }
 
522
                fprintf(matrix_file, "\n");
 
523
        }
 
524
        fclose(matrix_file);
 
525
        
 
526
        delete g_transMatrix;
 
527
}