-
Notifications
You must be signed in to change notification settings - Fork 0
/
asymptotic.c
120 lines (98 loc) · 4.18 KB
/
asymptotic.c
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
/* Copyright (c) 2014, Giuseppe Argentieri <[email protected]>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*
*/
#include "funcs.h"
#include "initial.h"
#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_ieee_utils.h>
int main ( int argc, char* argv[] )
{
gsl_ieee_env_setup () ; /* read GSL_IEEE_MODE */
printf("Temperature: %g\n\n", T*HBAR*Delta/BOLTZ) ;
double beta = 1.0/T ; /* Boltzmann factor: beta */
double omega_1 = gsl_hypot(OMEGA,D) ; /* omega' */
struct f_params params;
params.omega_c = omega_c ;
params.beta = beta ;
params.Omega = OMEGA ;
params.omega_1 = omega_1 ;
params.alpha = alpha ;
int status1 = save_integrals ( ¶ms ) ;
int status2 = save_matrices ( ¶ms ) ;
/* read the Redfield matrix from a file */
gsl_matrix* red_m = gsl_matrix_calloc ( 4, 4 ) ;
int status3 = mat_read ( red_m, "REDFIELD_MATRIX" ) ;
/* read the CP matrix from a file */
gsl_matrix* cp_m = gsl_matrix_calloc ( 4, 4 ) ;
int status4 = mat_read ( cp_m, "CP_MATRIX" ) ;
/* Hamiltonian generator */
const double om[] = { 0 , 0 , 0 , omega_1/2 } ;
gsl_matrix* H = gsl_matrix_calloc ( 4, 4 ) ;
int status5 = ham_gen ( H, om ) ;
/*
* Find the stationary state by solving the linear systems
* and the associated stationary currents
*/
gsl_vector* req_red = gsl_vector_calloc (4) ;
gsl_vector* req_cp = gsl_vector_calloc (4) ;
printf("REDFIELD DYNAMICS\n") ;
int status6 = stationary ( red_m , req_red ) ;
printf("Stationary state: ( %.1f , %.9f , %.9f , %.9f )\n",
VECTOR(req_red,0), VECTOR(req_red,1),
VECTOR(req_red,2), VECTOR(req_red,3) ) ;
printf("Stationary normalized (I/I0) DC current: %.9f\n\n",
-VECTOR(req_red,3)*OMEGA/omega_1) ;
printf("CP DYNAMICS\n") ;
int status7 = stationary ( cp_m , req_cp ) ;
printf("Stationary state: ( %.1f , %.9f , %.9f , %.9f )\n",
VECTOR(req_cp,0), VECTOR(req_cp,1),
VECTOR(req_cp,2), VECTOR(req_cp,3) ) ;
printf("Stationary normalized (I/I0) DC current: %.9f\n\n",
-VECTOR(req_cp,3)*OMEGA/omega_1) ;
/* Save the stationary states into files */
FILE* f = fopen ( "RED_STATIONARY.dat", "w" ) ;
gsl_vector_fprintf ( f, req_red, "%.9f" ) ;
FILE* g = fopen ( "CP_STATIONARY.dat", "w" ) ;
gsl_vector_fprintf ( g, req_cp, "%.9f" ) ;
fclose (f) ; fclose (g) ;
/* polarization */
double D30 = gsl_matrix_get(cp_m,3,0) ;
double D33 = gsl_matrix_get(cp_m,3,3) ;
double pol = -D30/D33 ;
printf("n-Polarization of the CP dynamics -D30/D33: %.9f\n", pol ) ;
/* polarization through the formula (for checking) */
double P = polarization ( ¶ms, omega_c ) ;
printf("\n") ;
printf("n-Polarization through the given formula: %.9f\n", P ) ;
printf("Stationary normalized (I/I0) DC current: %.9f\n", -P*OMEGA/omega_1 ) ;
/* free memory for matrices */
gsl_matrix_free(red_m) ;
gsl_matrix_free(cp_m) ;
return status1 + status2 + status3 + status4 + status5 + status6 + status7 ;
}