-
Notifications
You must be signed in to change notification settings - Fork 0
/
mkernel.wxm
128 lines (83 loc) · 3.25 KB
/
mkernel.wxm
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
/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 19.07.0 ] */
/* [wxMaxima: title start ]
Finds the memory kernel for the GBD simulations
[wxMaxima: title end ] */
/* [wxMaxima: comment start ]
Author: Andres Cordoba
Copyright (2022) Andres Cordoba
This script is distributed under the MIT License.
email: andcorduri@gmail.com
[wxMaxima: comment end ] */
/* [wxMaxima: section start ]
Functions for analysis
[wxMaxima: section end ] */
/* [wxMaxima: input start ] */
coeffs1(p,x):=block([l], l : [],
for i from 0 thru hipow(p,x)
do (l : cons(coeff(p,x,i),l)), l)$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
poly5(λ,H,ζ):=block([n],n: length(λ),
ζ*product(s+1/λ[j],j,1,n)+sum(H[j]*product(
if notequal(i,j) then
s+1/λ[i] else 1,i,1,n),j,1,n))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
poly4(n,ζ):=(μ*product(s+1/Λ[j],j,1,n)-sum(c[j]*product(
if notequal(i,j) then
s+1/Λ[i] else 1,i,1,n),j,1,n))*ζ$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
poly6(λ):=block([n],n:length(λ),product(s+1/λ[j],j,1,n))$
/* [wxMaxima: input end ] */
/* [wxMaxima: section start ]
System specification
[wxMaxima: section end ] */
/* [wxMaxima: input start ] */
input:read_list("input.csv", 'csv)$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
if input[4] > 0 then λin: makelist(input[i],i,6+input[4]+1,6+2*input[4]) else λin: []$
if input[4] > 0 then Hin: makelist(input[i],i,6,6+input[4]-1) else Hin: []$
ζin:input[6+2*input[4]+2]$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
write_data ([λin,Hin,ζin], "GsP.csv")$
/* [wxMaxima: input end ] */
/* [wxMaxima: section start ]
Calculate inertialess memory function
[wxMaxima: section end ] */
/* [wxMaxima: input start ] */
poly5eval:expand(poly5(λin,Hin,ζin))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
solΛ:allroots(poly5eval)$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
solΛ2:makelist(Λ[i]=at(-1/s,solΛ[i]),i,1,length(solΛ))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
poly4eval:expand(poly4(input[4],ζin))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
coeffs4:at(coeffs1(poly4eval,s),solΛ2)$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
poly6eval:expand(poly6(λin))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
coeffs6:coeffs1(poly6eval,s)$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
ecs:makelist(coeffs4[i]=coeffs6[i],i,1,length(coeffs6))$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
ratprint:false$
solc:solve(ecs),numer$
/* [wxMaxima: input end ] */
/* [wxMaxima: input start ] */
if input[4] > 0 then write_data ([at(makelist(Λ[i],i,1,length(solΛ)),solΛ2),at(makelist(c[i],i,1,length(solc[1])-1),solc[1])], "kernelp.csv") else write_data ([], "cs.csv")$
/* [wxMaxima: input end ] */
/* Old versions of Maxima abort on loading files that end in a comment. */
"Created with wxMaxima 19.07.0"$