-
Notifications
You must be signed in to change notification settings - Fork 0
/
ext_inputs.c
173 lines (153 loc) · 4.51 KB
/
ext_inputs.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
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
#include "ext_inputs.h"
struct pulse *new_pulse(int N, double t, double d, double q, double input,
double m)
{
struct pulse *p;
p = (struct pulse *) malloc(sizeof(struct pulse));
p->time_onset = t;
p->duration = d;
p->theta = q; /* q is given in degrees! */
p->Imax = input;
p->concentration = m;
p->switched_on = 0;
p->precomputed_value = precompute_vector_pulses(N, p);
p->next = NULL;
return p;
}
struct pulse_2D *new_pulse_2D(int N1, int N2, double t, double d, double q1,
double q2, double input, double m)
{
struct pulse_2D *p;
p = (struct pulse_2D *) malloc(sizeof(struct pulse_2D));
p->time_onset = t;
p->duration = d;
p->theta1 = q1;
p->theta2 = q2;
p->Imax = input;
p->concentration = m;
p->switched_on = 0;
p->precomputed_value = precompute_vector_pulses_2D(N1, N2, p);
p->next = NULL;
return p;
}
double *precompute_vector_pulses(int N, struct pulse *plse)
{
double *v = malloc(N * sizeof(double));
double angle[N];
double theta_rad = plse->theta * M_PI / 180.;
for (int j = 0; j < N; j++) {
/* TODO avoid computing angle each time a new pulse is added */
angle[j] = M_PI * (-1. / 2. + (double) j / (double) N);
v[j] = plse->Imax * exp(plse->concentration
* (cos(2.0 * (angle[j] - theta_rad)) - 1));
}
return v;
}
double *precompute_vector_pulses_2D(int N1, int N2, struct pulse_2D *plse)
{
double *v = malloc(N1 * N2 * sizeof(double));
double angle1[N1];
double angle2[N2];
double theta1_rad = plse->theta1 * M_PI / 180.;
double theta2_rad = plse->theta2 * M_PI / 180.;
double dist;
/* TODO avoid computing angle each time a new pulse is added */
for (int i = 0; i < N1; i++)
angle1[i] = M_PI * (-1. / 2. + (double) i / (double) N1);
for (int i = 0; i < N2; i++)
angle2[i] = M_PI * (-1. / 2. + (double) i / (double) N2);
for (int i = 0; i < N1; i++) {
for (int j = 0; j < N2; j++) {
dist =
distance_on_torus(angle1[i], angle2[j], theta1_rad, theta2_rad);
v[N2 * i + j] =
plse->Imax * exp(plse->concentration * (cos(2.0 * dist) - 1));
}
}
return v;
}
struct pulse *append_pulse(struct pulse *plse, struct pulse *newplse)
{
struct pulse *p;
/* this is the first item in the list */
if (plse == NULL) {
return newplse;
} else {
/* walk to the end of the list */
for (p = plse; p->next != NULL; p = p->next)
/* do nothing */ ;
p->next = newplse;
return plse;
}
}
struct pulse_2D *append_pulse_2D(struct pulse_2D *plse,
struct pulse_2D *newplse)
{
struct pulse_2D *p;
/* this is the first item in the list */
if (plse == NULL) {
return newplse;
} else {
/* walk to the end of the list */
for (p = plse; p->next != NULL; p = p->next)
/* do nothing */ ;
p->next = newplse;
return plse;
}
}
void print_pulses(struct pulse *plse)
{
printf("\nExternal inputs...\n");
for (; plse != NULL; plse = plse->next) {
printf("% 8.2f % 8.2f % 8.2f % 8.2f % 8.2f\n",
plse->time_onset, plse->duration, plse->theta,
plse->Imax, plse->concentration);
}
}
void print_pulses_2D(struct pulse_2D *plse)
{
printf("\nExternal inputs...\n");
for (; plse != NULL; plse = plse->next) {
printf("% 8.2f % 8.2f % 8.2f % 8.2f % 8.2f % 8.2f\n",
plse->time_onset, plse->duration,
plse->theta1, plse->theta2, plse->Imax, plse->concentration);
}
}
double max_amplitude_pulse(struct pulse *plse)
{
double tmp = 0;
for (; plse != NULL; plse = plse->next) {
if (plse->Imax > tmp) {
tmp = plse->Imax;
}
}
return tmp;
}
double max_amplitude_pulse_2D(struct pulse_2D *plse)
{
double tmp = 0;
for (; plse != NULL; plse = plse->next) {
if (plse->Imax > tmp) {
tmp = plse->Imax;
}
}
return tmp;
}
void free_extinputs(struct pulse *plse)
{
struct pulse *next;
for (; plse != NULL; plse = next) {
next = plse->next;
free(plse->precomputed_value);
free(plse);
}
}
void free_extinputs_2D(struct pulse_2D *plse)
{
struct pulse_2D *next;
for (; plse != NULL; plse = next) {
next = plse->next;
free(plse->precomputed_value);
free(plse);
}
}