software/simplesolv
1
Fork 0
Code Issues Pull Requests Releases Activity

Update to v0.4

Browse Source 
feature: compute the 2-point correlation function in easyeq.

  feature: compute the Fourier transform of the 2-point correlation function
           in anyeq and easyeq.

  feature: compute the local maximum of the 2-point correlation function and
           its Fourier transform.

  feature: compute the compressibility for anyeq.

  feature: allow for linear spacing of rho's.

  feature: print the scattering length.

  change: ux and uk now return real numbers.

  fix: error in the computation of the momentum distribution: wrong
       definition of delta functions.

  fix: various minor bugs.

  optimization: assign explicit types to variables.
This commit is contained in:
Ian Jauslin 2023-02-26 18:36:05 -05:00
parent e72af82c3e
commit c1b477a1b2
26 changed files with 4310 additions and 841 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
Changelog Normal file
View File

@ -0,0 +1,25 @@
0.4:
* feature: compute the 2-point correlation function in easyeq.
* feature: compute the Fourier transform of the 2-point correlation function
in anyeq and easyeq.
* feature: compute the local maximum of the 2-point correlation function and
its Fourier transform.
* feature: compute the compressibility for anyeq.
* feature: allow for linear spacing of rho's.
* feature: print the scattering length.
* change: ux and uk now return real numbers.
* fix: error in the computation of the momentum distribution: wrong
definition of delta functions.
* fix: various minor bugs.
* optimization: assign explicit types to variables.

2
NOTICE
View File

@ -1,4 +1,4 @@
simplesolv: simplesolv:
Licensed under the Apache 2.0 License (see LICENSE for details) Licensed under the Apache 2.0 License (see LICENSE for details)
Copyright Ian Jauslin 2021 Copyright Ian Jauslin 2021-2023

11
README
View File

@ -29,7 +29,9 @@ published in the following papers:
* E.A. Carlen, I. Jauslin, E.H. Lieb - Analysis of a simple equation for the * E.A. Carlen, I. Jauslin, E.H. Lieb - Analysis of a simple equation for the
ground state energy of the Bose gas II: Monotonicity, Convexity and ground state energy of the Bose gas II: Monotonicity, Convexity and
Condensate Fraction, to appear in the SIAM Journal on Mathematical Analysis Condensate Fraction, SIAM Journal on Mathematical Analysis, volume 53,
number 5, pages 5322-5360, 2021
https://doi.org/10.1137/20M13768
https://arxiv.org/abs/2010.13882 https://arxiv.org/abs/2010.13882
http://ian.jauslin.org/publications/20cjl http://ian.jauslin.org/publications/20cjl
@ -40,6 +42,13 @@ published in the following papers:
https://arxiv.org/abs/2011.10869 https://arxiv.org/abs/2011.10869
http://ian.jauslin.org/publications/20chjl http://ian.jauslin.org/publications/20chjl
* I. Jauslin - The Simplified approach to the Bose gas without translation
invariance, 2023
http://ian.jauslin.org/publications/23j
* I. Jauslin - Evidence of a liquid phase in repulsive Bosonic systems using
the Simplified approach to the Bose gas, 2023
http://ian.jauslin.org/publications/23j_b
Dependencies: Dependencies:

2
doc/Makefile
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.

7
doc/bibliography/bibliography.tex
View File

@ -4,10 +4,13 @@ doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRevA.103.053309}{10.1103
\bibitem[CJL20]{CJL20}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a simple equation for the ground state energy of the Bose gas}, Pure and Applied Analysis, volume~\-2, issue~\-3, pages~\-659-684, 2020,\par\penalty10000 \bibitem[CJL20]{CJL20}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a simple equation for the ground state energy of the Bose gas}, Pure and Applied Analysis, volume~\-2, issue~\-3, pages~\-659-684, 2020,\par\penalty10000
doi:{\tt\color{blue}\href{http://dx.doi.org/10.2140/paa.2020.2.659}{10.2140/paa.2020.2.659}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1912.04987}{1912.04987}}.\par\medskip doi:{\tt\color{blue}\href{http://dx.doi.org/10.2140/paa.2020.2.659}{10.2140/paa.2020.2.659}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/1912.04987}{1912.04987}}.\par\medskip
\bibitem[CJL20b]{CJL20b}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a simple equation for the ground state of the Bose gas II: Monotonicity, Convexity and Condensate Fraction}, 2020, to appear in the SIAM journal of Mathematical Analysis,\par\penalty10000 \bibitem[CJL21]{CJL21}E.A. Carlen, I. Jauslin, E.H. Lieb - {\it Analysis of a Simple Equation for the Ground State of the Bose Gas II: Monotonicity, Convexity, and Condensate Fraction}, SIAM Journal on Mathematical Analysis, volume~\-53, number~\-5, pages~\-5322-5360, 2021,\par\penalty10000
arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/2010.13882}{2010.13882}}.\par\medskip doi:{\tt\color{blue}\href{http://dx.doi.org/10.1137/20M1376820}{10.1137/20M1376820}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/2010.13882}{2010.13882}}.\par\medskip
\bibitem[DLMF]{DLMF1.1.3}F.W.J. Olver, A.B. Olde Daalhuis, D.W. Lozier, B.I. Schneider, R.F. Boisvert, C.W. Clark, B.R. Miller, B.V. Saunders, H.S. Cohl, M.A. McClain (editors) - {\it NIST Digital Library of Mathematical Functions}, Release~\-1.1.3 of~\-2021-09-15, 2021.\par\medskip \bibitem[DLMF]{DLMF1.1.3}F.W.J. Olver, A.B. Olde Daalhuis, D.W. Lozier, B.I. Schneider, R.F. Boisvert, C.W. Clark, B.R. Miller, B.V. Saunders, H.S. Cohl, M.A. McClain (editors) - {\it NIST Digital Library of Mathematical Functions}, Release~\-1.1.3 of~\-2021-09-15, 2021.\par\medskip
\bibitem[Ja23]{Ja23}I. Jauslin - {\it The Simplified approach to the Bose gas without translation invariance}, 2023\par\penalty10000
{\tt\color{blue}\href{http://ian.jauslin.org/publications/23j/}{http://ian.jauslin.org/publications/23j/}}.\par\medskip
\bibitem[Ta87]{Ta87}Y. Taguchi - {\it Fourier coefficients of periodic functions of Gevrey classes and ultradistributions}, Yokohama Mathematical Journal, volume~\-35, pages~\-51-60, 1987.\par\medskip \bibitem[Ta87]{Ta87}Y. Taguchi - {\it Fourier coefficients of periodic functions of Gevrey classes and ultradistributions}, Yokohama Mathematical Journal, volume~\-35, pages~\-51-60, 1987.\par\medskip

2
doc/libs/code.sty
View File

@ -1,4 +1,4 @@
%% Copyright 2021 Ian Jauslin %% Copyright 2021-2023 Ian Jauslin
%% %%
%% Licensed under the Apache License, Version 2.0 (the "License"); %% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License. %% you may not use this file except in compliance with the License.

2
doc/libs/dlmf.sty
View File

@ -1,4 +1,4 @@
%% Copyright 2021 Ian Jauslin %% Copyright 2021-2023 Ian Jauslin
%% %%
%% Licensed under the Apache License, Version 2.0 (the "License"); %% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License. %% you may not use this file except in compliance with the License.

53
doc/libs/ian.cls
View File

@ -219,6 +219,7 @@
% counters % counters
\stepcounter{subsectioncount} \stepcounter{subsectioncount}
\setcounter{subsubsectioncount}{0} \setcounter{subsubsectioncount}{0}
\setcounter{subsubsubsectioncount}{0}
\ifsubsectionsineq \ifsubsectionsineq
\setcounter{seqcount}0 \setcounter{seqcount}0
\fi \fi
@ -268,6 +269,7 @@
\def\subsubsection#1{ \def\subsubsection#1{
% counters % counters
\stepcounter{subsubsectioncount} \stepcounter{subsubsectioncount}
\setcounter{subsubsubsectioncount}{0}
% space before subsubsection (if not first) % space before subsubsection (if not first)
\ifnum\thesubsubsectioncount>1 \ifnum\thesubsubsectioncount>1
@ -302,6 +304,49 @@
\medskip\penalty10000 \medskip\penalty10000
} }
%% subsubsubsection
% counter
\newcounter{subsubsubsectioncount}
% space before subsubsubsection
\newlength\subsubsubsecskip
\setlength\subsubsubsecskip{20pt}
\def\subsubsubsection#1{
% counters
\stepcounter{subsubsubsectioncount}
% space before subsubsubsection (if not first)
\ifnum\thesubsubsubsectioncount>1
\vskip\subsubsubsecskip
\penalty-500
\fi
% hyperref anchor
\hrefanchor
% define tag (for \label)
\xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\thesubsubsubsectioncount}
% get widths
\def\@secnum{{\bf\hskip1.5cm\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\thesubsubsubsectioncount.\hskip5pt}}
\settowidth\secnumwidth{\@secnum}
\setlength\sectitlewidth\textwidth
\addtolength\sectitlewidth{-\secnumwidth}
% print name
\parbox{\textwidth}{
\@secnum
\parbox[t]{\sectitlewidth}{\bf #1}}
% write to table of contents
\iftoc
% save lncount in aux variable which is written to toc
\immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@subsubsubsec.\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\thesubsubsubsectioncount\endcsname{\thelncount}}
\write\tocoutput{\noexpand\tocsubsubsubsection{#1}{\thepage}}
\fi
% space
\par\penalty10000
\medskip\penalty10000
}
%% itemize %% itemize
\newlength\itemizeskip \newlength\itemizeskip
% left margin for items % left margin for items
@ -638,6 +683,7 @@
\stepcounter{tocsectioncount} \stepcounter{tocsectioncount}
\setcounter{tocsubsectioncount}{0} \setcounter{tocsubsectioncount}{0}
\setcounter{tocsubsubsectioncount}{0} \setcounter{tocsubsubsectioncount}{0}
\setcounter{tocsubsubsubsectioncount}{0}
% write % write
\smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par \smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par
} }
@ -645,6 +691,7 @@
\def\tocsubsection #1#2{ \def\tocsubsection #1#2{
\stepcounter{tocsubsectioncount} \stepcounter{tocsubsectioncount}
\setcounter{tocsubsubsectioncount}{0} \setcounter{tocsubsubsectioncount}{0}
\setcounter{tocsubsubsubsectioncount}{0}
% write % write
{\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par {\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
} }
@ -654,6 +701,12 @@
% write % write
{\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par {\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
} }
\newcounter{tocsubsubsubsectioncount}
\def\tocsubsubsubsection #1#2{
\stepcounter{tocsubsubsubsectioncount}
% write
{\hskip30pt\hyperlink{ln.\csname toc@subsubsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount.\thetocsubsubsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount.\thetocsubsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
}
\def\tocappendices{ \def\tocappendices{
\medskip \medskip
\setcounter{tocsectioncount}0 \setcounter{tocsectioncount}0

2
doc/libs/iantheo.sty
View File

@ -1,4 +1,4 @@
%% Copyright 2021 Ian Jauslin %% Copyright 2021-2023 Ian Jauslin
%% %%
%% Licensed under the Apache License, Version 2.0 (the "License"); %% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License. %% you may not use this file except in compliance with the License.

2
doc/libs/largearray.sty
View File

@ -1,4 +1,4 @@
%% Copyright 2021 Ian Jauslin %% Copyright 2021-2023 Ian Jauslin
%% %%
%% Licensed under the Apache License, Version 2.0 (the "License"); %% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License. %% you may not use this file except in compliance with the License.

12
doc/libs/point.sty
View File

@ -1,4 +1,4 @@
%% Copyright 2021 Ian Jauslin %% Copyright 2021-2023 Ian Jauslin
%% %%
%% Licensed under the Apache License, Version 2.0 (the "License"); %% Licensed under the Apache License, Version 2.0 (the "License");
%% you may not use this file except in compliance with the License. %% you may not use this file except in compliance with the License.
@ -33,12 +33,15 @@
\newif\ifresetatsubsubsection \newif\ifresetatsubsubsection
\DeclareOption{reset_at_subsubsection}{\resetatsubsubsectiontrue} \DeclareOption{reset_at_subsubsection}{\resetatsubsubsectiontrue}
\DeclareOption{no_reset_at_subsubsection}{\resetatsubsubsectionfalse} \DeclareOption{no_reset_at_subsubsection}{\resetatsubsubsectionfalse}
\newif\ifresetatsubsubsubsection
\DeclareOption{reset_at_subsubsubsection}{\resetatsubsubsubsectiontrue}
\DeclareOption{no_reset_at_subsubsubsection}{\resetatsubsubsubsectionfalse}
\newif\ifresetattheo \newif\ifresetattheo
\DeclareOption{reset_at_theo}{\resetattheotrue} \DeclareOption{reset_at_theo}{\resetattheotrue}
\DeclareOption{no_reset_at_theo}{\resetattheofalse} \DeclareOption{no_reset_at_theo}{\resetattheofalse}
\def\point@defaultoptions{ \def\point@defaultoptions{
\ExecuteOptions{reset_at_section, reset_at_subsection, reset_at_subsubsection, no_reset_at_theo} \ExecuteOptions{reset_at_section, reset_at_subsection, reset_at_subsubsection, reset_at_subsubsubsection, no_reset_at_theo}
\ProcessOptions \ProcessOptions
%% reset at every new section %% reset at every new section
@ -56,6 +59,11 @@
\let\point@oldsubsubsection\subsubsection \let\point@oldsubsubsection\subsubsection
\gdef\subsubsection{\resetpointcounter\point@oldsubsubsection} \gdef\subsubsection{\resetpointcounter\point@oldsubsubsection}
\fi \fi
%% reset at every new subsubsubsection
\ifresetatsubsubsubsection
\let\point@oldsubsubsubsection\subsubsubsection
\gdef\subsubsubsection{\resetpointcounter\point@oldsubsubsubsection}
\fi
%% reset at every new theorem %% reset at every new theorem
\ifresetattheo \ifresetattheo

1110
doc/simplesolv-doc.tex
View File

File diff suppressed because it is too large Load Diff

1505
src/anyeq.jl
View File

File diff suppressed because it is too large Load Diff

323
src/chebyshev.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,7 +13,15 @@
## limitations under the License. ## limitations under the License.
# Chebyshev expansion # Chebyshev expansion
@everywhere function chebyshev(a,taus,weights,P,N,J,nu) @everywhere function chebyshev(
a::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=zeros(Float64,J*(P+1)) out=zeros(Float64,J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
@ -27,7 +35,15 @@
end end
# evaluate function from Chebyshev expansion # evaluate function from Chebyshev expansion
@everywhere function chebyshev_eval(Fa,x,taus,chebyshev,P,J,nu) @everywhere function chebyshev_eval(
Fa::Array{Float64,1},
x::Float64,
taus::Array{Float64,1},
chebyshev::Array{Polynomial,1},
P::Int64,
J::Int64,
nu::Int64
)
# change variable # change variable
tau=(1-x)/(1+x) tau=(1-x)/(1+x)
@ -48,7 +64,11 @@ end
# convolution # convolution
# input the Chebyshev expansion of a and b, as well as the A matrix # input the Chebyshev expansion of a and b, as well as the A matrix
@everywhere function conv_chebyshev(Fa,Fb,A) @everywhere function conv_chebyshev(
Fa::Array{Float64,1},
Fb::Array{Float64,1},
A::Array{Array{Float64,2},1}
)
out=zeros(Float64,length(A)) out=zeros(Float64,length(A))
for i in 1:length(A) for i in 1:length(A)
out[i]=dot(Fa,A[i]*Fb) out[i]=dot(Fa,A[i]*Fb)
@ -57,12 +77,27 @@ end
end end
# <ab> # <ab>
@everywhere function avg_chebyshev(Fa,Fb,A0) @everywhere function avg_chebyshev(
Fa::Array{Float64,1},
Fb::Array{Float64,1},
A0::Float64
)
return dot(Fa,A0*Fb) return dot(Fa,A0*Fb)
end end
# 1_n * a # 1_n * a
@everywhere function conv_one_chebyshev(n,zetapp,Fa,A,taus,weights,P,N,J,nu1) @everywhere function conv_one_chebyshev(
n::Int64,
zetapp::Int64,
Fa::Array{Float64,1},
A::Array{Array{Float64,2},1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu1::Int64
)
out=zeros(Float64,N*J) out=zeros(Float64,N*J)
for m in 1:N*J for m in 1:N*J
for l in 0:P for l in 0:P
@ -74,7 +109,15 @@ end
return out return out
end end
# a * v # a * v
@everywhere function conv_v_chebyshev(a,Upsilon,k,taus,weights,N,J) @everywhere function conv_v_chebyshev(
a::Array{Float64,1},
Upsilon::Array{Array{Float64,1},1},
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=zeros(Float64,J*N) out=zeros(Float64,J*N)
for i in 1:J*N for i in 1:J*N
for zetap in 0:J-1 for zetap in 0:J-1
@ -86,7 +129,16 @@ end
end end
return out return out
end end
@everywhere function conv_one_v_chebyshev(n,zetap,Upsilon,k,taus,weights,N,J) @everywhere function conv_one_v_chebyshev(
n::Int64,
zetap::Int64,
Upsilon::Array{Array{Float64,1},1},
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=zeros(Float64,J*N) out=zeros(Float64,J*N)
xj=weights[1][n] xj=weights[1][n]
for i in 1:J*N for i in 1:J*N
@ -96,7 +148,14 @@ end
end end
# <av> # <av>
@everywhere function avg_v_chebyshev(a,Upsilon0,k,taus,weights,N,J) @everywhere function avg_v_chebyshev(a,
Upsilon0::Array{Float64,1},
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=0. out=0.
for zetap in 0:J-1 for zetap in 0:J-1
for j in 1:N for j in 1:N
@ -107,13 +166,28 @@ end
return out return out
end end
# <1_nv> # <1_nv>
@everywhere function avg_one_v_chebyshev(n,zetap,Upsilon0,k,taus,weights,N) @everywhere function avg_one_v_chebyshev(
n::Int64,
zetap::Int64,
Upsilon0::Array{Float64,1},
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64
)
xj=weights[1][n] xj=weights[1][n]
return (taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon0[zetap*N+n] return (taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon0[zetap*N+n]
end end
# compute \int dq dxi u1(k-xi)u2(q)u3(xi)u4(k-q)u5(xi-q) # compute \int dq dxi u1(k-xi)u2(q)u3(xi)u4(k-q)u5(xi-q)
@everywhere function double_conv_S_chebyshev(FU1,FU2,FU3,FU4,FU5,Abar) @everywhere function double_conv_S_chebyshev(
FU1::Array{Float64,1},
FU2::Array{Float64,1},
FU3::Array{Float64,1},
FU4::Array{Float64,1},
FU5::Array{Float64,1},
Abar::Array{Float64,5}
)
out=zeros(Float64,length(Abar)) out=zeros(Float64,length(Abar))
for i in 1:length(Abar) for i in 1:length(Abar)
for j1 in 1:length(FU1) for j1 in 1:length(FU1)
@ -133,7 +207,17 @@ end
# compute A # compute A
@everywhere function Amat(k,weights,taus,T,P,N,J,nua,nub) @everywhere function Amat(
k::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
taus::Array{Float64,1},
T::Array{Polynomial,1},
P::Int64,
N::Int64,
J::Int64,
nua::Int64,
nub::Int64
)
out=Array{Array{Float64,2},1}(undef,J*N) out=Array{Array{Float64,2},1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=zeros(Float64,J*(P+1),J*(P+1)) out[i]=zeros(Float64,J*(P+1),J*(P+1))
@ -152,7 +236,11 @@ end
end end
# compute Upsilon # compute Upsilon
@everywhere function Upsilonmat(k,v,weights) @everywhere function Upsilonmat(
k::Array{Float64,1},
v::Function,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
out=Array{Array{Float64,1},1}(undef,length(k)) out=Array{Array{Float64,1},1}(undef,length(k))
for i in 1:length(k) for i in 1:length(k)
out[i]=Array{Float64,1}(undef,length(k)) out[i]=Array{Float64,1}(undef,length(k))
@ -162,7 +250,11 @@ end
end end
return out return out
end end
@everywhere function Upsilon0mat(k,v,weights) @everywhere function Upsilon0mat(
k::Array{Float64,1},
v::Function,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
out=Array{Float64,1}(undef,length(k)) out=Array{Float64,1}(undef,length(k))
for j in 1:length(k) for j in 1:length(k)
out[j]=2*k[j]*v(k[j]) out[j]=2*k[j]*v(k[j])
@ -171,17 +263,36 @@ end
end end
# alpha_- # alpha_-
@everywhere function alpham(k,t) @everywhere function alpham(
k::Float64,
t::Float64
)
return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t)) return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t))
end end
# alpha_+ # alpha_+
@everywhere function alphap(k,t) @everywhere function alphap(
k::Float64,
t::Float64
)
return (1-abs(k-(1-t)/(1+t)))/(1+abs(k-(1-t)/(1+t))) return (1-abs(k-(1-t)/(1+t)))/(1+abs(k-(1-t)/(1+t)))
end end
# compute \bar A # compute \bar A
@everywhere function barAmat(k,weights,taus,T,P,N,J,nu1,nu2,nu3,nu4,nu5) @everywhere function barAmat(
k::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}},
taus::Array{Float64,1},
T::Array{Polynomial,1},
P::Int64,
N::Int64,
J::Int64,
nu1::Int64,
nu2::Int64,
nu3::Int64,
nu4::Int64,
nu5::Int64
)
out=zeros(Float64,J*(P+1),J*(P+1),J*(P+1),J*(P+1),J*(P+1)) out=zeros(Float64,J*(P+1),J*(P+1),J*(P+1),J*(P+1),J*(P+1))
for zeta1 in 0:J-1 for zeta1 in 0:J-1
for n1 in 0:P for n1 in 0:P
@ -211,27 +322,107 @@ end
return out return out
end end
@everywhere function barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5) @everywhere function barAmat_int1(tau,
k::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu1::Int64,
nu2::Int64,
nu3::Int64,
nu4::Int64,
nu5::Int64,
zeta1::Int64,
zeta2::Int64,
zeta3::Int64,
zeta4::Int64,
zeta5::Int64,
n1::Int64,
n2::Int64,
n3::Int64,
n4::Int64,
n5::Int64
)
if(alpham(k,tau)<taus[zeta2+2] && alphap(k,tau)>taus[zeta2+1]) if(alpham(k,tau)<taus[zeta2+2] && alphap(k,tau)>taus[zeta2+1])
return 2*(1-tau)/(1+tau)^(3-nu1)*T[n1+1]((2*tau-(taus[zeta1+1]+taus[zeta1+2]))/(taus[zeta1+2]-taus[zeta1+1]))*integrate_legendre(sigma->barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5),max(taus[zeta2+1],alpham(k,tau)),min(taus[zeta2+2],alphap(k,tau)),weights) return 2*(1-tau)/(1+tau)^(3-nu1)*T[n1+1]((2*tau-(taus[zeta1+1]+taus[zeta1+2]))/(taus[zeta1+2]-taus[zeta1+1]))*integrate_legendre(sigma->barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5),max(taus[zeta2+1],alpham(k,tau)),min(taus[zeta2+2],alphap(k,tau)),weights)
else else
return 0. return 0.
end end
end end
@everywhere function barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5) @everywhere function barAmat_int2(tau,
sigma::Float64,
k::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu2::Int64,
nu3::Int64,
nu4::Int64,
nu5::Int64,
zeta2::Int64,
zeta3::Int64,
zeta4::Int64,
zeta5::Int64,
n2::Int64,
n3::Int64,
n4::Int64,
n5::Int64
)
return 2*(1-sigma)/(1+sigma)^(3-nu2)*T[n2+1]((2*sigma-(taus[zeta2+1]+taus[zeta2+2]))/(taus[zeta2+2]-taus[zeta2+1]))*integrate_legendre(taup->barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5),taus[zeta3+1],taus[zeta3+2],weights) return 2*(1-sigma)/(1+sigma)^(3-nu2)*T[n2+1]((2*sigma-(taus[zeta2+1]+taus[zeta2+2]))/(taus[zeta2+2]-taus[zeta2+1]))*integrate_legendre(taup->barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5),taus[zeta3+1],taus[zeta3+2],weights)
end end
@everywhere function barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5) @everywhere function barAmat_int3(tau,
sigma::Float64,
taup::Float64,
k::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu3::Int64,
nu4::Int64,
nu5::Int64,
zeta3::Int64,
zeta4::Int64,
zeta5::Int64,
n3::Int64,
n4::Int64,
n5::Int64
)
if(alpham(k,taup)<taus[zeta4+2] && alphap(k,taup)>taus[zeta4+1]) if(alpham(k,taup)<taus[zeta4+2] && alphap(k,taup)>taus[zeta4+1])
return 2*(1-taup)/(1+taup)^(3-nu3)*T[n3+1]((2*taup-(taus[zeta3+1]+taus[zeta3+2]))/(taus[zeta3+2]-taus[zeta3+1]))*integrate_legendre(sigmap->barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5),max(taus[zeta4+1],alpham(k,taup)),min(taus[zeta4+2],alphap(k,taup)),weights) return 2*(1-taup)/(1+taup)^(3-nu3)*T[n3+1]((2*taup-(taus[zeta3+1]+taus[zeta3+2]))/(taus[zeta3+2]-taus[zeta3+1]))*integrate_legendre(sigmap->barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5),max(taus[zeta4+1],alpham(k,taup)),min(taus[zeta4+2],alphap(k,taup)),weights)
else else
return 0. return 0.
end end
end end
@everywhere function barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5) @everywhere function barAmat_int4(tau,
sigma::Float64,
taup::Float64,
sigmap::Float64,
k::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu4::Int64,
nu5::Int64,
zeta4::Int64,
zeta5::Int64,
n4::Int64,
n5::Int64
)
return 2*(1-sigmap)/(1+sigmap)^(3-nu4)*T[n4+1]((2*sigma-(taus[zeta4+1]+taus[zeta4+2]))/(taus[zeta4+2]-taus[zeta4+1]))*integrate_legendre(theta->barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5),0.,2*pi,weights) return 2*(1-sigmap)/(1+sigmap)^(3-nu4)*T[n4+1]((2*sigma-(taus[zeta4+1]+taus[zeta4+2]))/(taus[zeta4+2]-taus[zeta4+1]))*integrate_legendre(theta->barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5),0.,2*pi,weights)
end end
@everywhere function barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5) @everywhere function barAmat_int5(tau,
sigma::Float64,
taup::Float64,
sigmap::Float64,
theta::Float64,
k::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu5::Int64,
zeta5::Int64,
n5::Int64
)
R=barAmat_R((1-sigma)/(1+sigma),(1-tau)/(1+tau),(1-sigmap)/(1+sigmap),(1-taup)/(1+taup),theta,k) R=barAmat_R((1-sigma)/(1+sigma),(1-tau)/(1+tau),(1-sigmap)/(1+sigmap),(1-taup)/(1+taup),theta,k)
if((1-R)/(1+R)<taus[zeta5+2] && (1-R)/(1+R)>taus[zeta5+1]) if((1-R)/(1+R)<taus[zeta5+2] && (1-R)/(1+R)>taus[zeta5+1])
return (2/(2+R))^nu5*T[n5+1]((2*(1-R)/(1+R)-(taus[zeta5+1]+taus[zeta5+2]))/(taus[zeta5+2]-taus[zeta5+1])) return (2/(2+R))^nu5*T[n5+1]((2*(1-R)/(1+R)-(taus[zeta5+1]+taus[zeta5+2]))/(taus[zeta5+2]-taus[zeta5+1]))
@ -240,13 +431,22 @@ end
end end
end end
# R(s,t,s',t,theta,k) # R(s,t,s',t,theta,k)
@everywhere function barAmat_R(s,t,sp,tp,theta,k) @everywhere function barAmat_R(
return sqrt(k^2*(s^2+t^2+sp^2+tp^2)-k^4-(s^2-t^2)*(sp^2-tp^2)-sqrt((4*k^2*s^2-(k^2+s^2-t^2)^2)*(4*k^2*sp^2-(k^2+sp^2-tp^2)^2))*cos(theta))/(sqrt(2)*k) s::Float64,
t::Float64,
sp::Float64,
tp::Float64,
theta::Float64,
k::Float64
)
return sqrt(k^2*(s^2+t^2+sp^2+tp^2)-k^4-(s^2-t^2)*(sp^2-tp^2)-sqrt((4*k^2*s^2-(k^2+s^2-t^2)^2)*(4*k^2*sp^2-(k^2+sp^2-tp^2)^2))*cos(theta))/(sqrt(2.)*k)
end end
# compute Chebyshev polynomials # compute Chebyshev polynomials
@everywhere function chebyshev_polynomials(P) @everywhere function chebyshev_polynomials(
T=Array{Polynomial}(undef,P+1) P::Int64
)
T=Array{Polynomial,1}(undef,P+1)
T[1]=Polynomial([1]) T[1]=Polynomial([1])
T[2]=Polynomial([0,1]) T[2]=Polynomial([0,1])
for n in 1:P-1 for n in 1:P-1
@ -258,7 +458,15 @@ end
end end
# compute \int f*u dk/(2*pi)^3 # compute \int f*u dk/(2*pi)^3
@everywhere function integrate_f_chebyshev(f,u,k,taus,weights,N,J) @everywhere function integrate_f_chebyshev(
f::Function,
u::Array{Float64,1},
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=0. out=0.
for zeta in 0:J-1 for zeta in 0:J-1
for i in 1:N for i in 1:N
@ -268,7 +476,15 @@ end
return out return out
end end
@everywhere function inverse_fourier_chebyshev(u,x,k,taus,weights,N,J) @everywhere function inverse_fourier_chebyshev(
u::Array{Float64,1},
x::Float64,
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=0. out=0.
for zeta in 0:J-1 for zeta in 0:J-1
for j in 1:N for j in 1:N
@ -277,3 +493,54 @@ end
end end
return out return out
end end
# compute B (for the computation of the fourier transform of the two-point correlation)
@everywhere function Bmat(
q::Float64,
k::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
taus::Array{Float64,1},
T::Array{Polynomial,1},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,1},1}(undef,J*N)
for i in 1:J*N
out[i]=zeros(Float64,J*(P+1))
for zeta in 0:J-1
for n in 0:P
out[i][zeta*(P+1)+n+1]=1/(8*pi^3*k[i]*q)*(betam(k[i],q)>taus[zeta+2] || betap(k[i],q)<taus[zeta+1] ? 0. : integrate_legendre(sigma->(1-sigma)/(1+sigma)^(3-nu)*T[n+1]((2*sigma-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],betam(k[i],q)),min(taus[zeta+2],betap(k[i],q)),weights))
end
end
end
return out
end
# beta_-
@everywhere function betam(
k::Float64,
q::Float64
)
return (1-k-q)/(1+k+q)
end
# beta_+
@everywhere function betap(
k::Float64,
q::Float64
)
return (1-abs(k-q))/(1+abs(k-q))
end
# mathfrak S (for the computation of the fourier transform of the two-point correlation)
@everywhere function chebyshev_frakS(
Ff::Array{Float64,1},
B::Array{Array{Float64,1},1}
)
out=zeros(Float64,length(B))
for i in 1:length(B)
out[i]=dot(Ff,B[i])
end
return out
end

1044
src/easyeq.jl
View File

File diff suppressed because it is too large Load Diff

50
src/integration.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,7 +13,12 @@
## limitations under the License. ## limitations under the License.
# approximate \int_a^b f using Gauss-Legendre quadratures # approximate \int_a^b f using Gauss-Legendre quadratures
@everywhere function integrate_legendre(f,a,b,weights) @everywhere function integrate_legendre(
f::Function,
a::Float64,
b::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
out=0 out=0
for i in 1:length(weights[1]) for i in 1:length(weights[1])
out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2) out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)
@ -21,7 +26,13 @@
return out return out
end end
# \int f*g where g is sampled at the Legendre nodes # \int f*g where g is sampled at the Legendre nodes
@everywhere function integrate_legendre_sampled(f,g,a,b,weights) @everywhere function integrate_legendre_sampled(
f::Function,
g::Array{Float64,1},
a::Float64,
b::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
out=0 out=0
for i in 1:length(weights[1]) for i in 1:length(weights[1])
out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i] out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i]
@ -31,7 +42,12 @@ end
# approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures # approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures
@everywhere function integrate_chebyshev(f,a,b,N) @everywhere function integrate_chebyshev(
f::Function,
a::Float64,
b::Float64,
N::Int64
)
out=0 out=0
for i in 1:N for i in 1:N
out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2) out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2)
@ -40,7 +56,11 @@ end
end end
# approximate \int_0^\infty dr f(r)*exp(-a*r) using Gauss-Chebyshev quadratures # approximate \int_0^\infty dr f(r)*exp(-a*r) using Gauss-Chebyshev quadratures
@everywhere function integrate_laguerre(f,a,weights_gL) @everywhere function integrate_laguerre(
f::Function,
a::Float64,
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}}
)
out=0. out=0.
for i in 1:length(weights_gL[1]) for i in 1:length(weights_gL[1])
out+=1/a*f(weights_gL[1][i]/a)*weights_gL[2][i] out+=1/a*f(weights_gL[1][i]/a)*weights_gL[2][i]
@ -49,10 +69,28 @@ end
end end
# Hann window # Hann window
@everywhere function hann(x,L) @everywhere function hann(
x::Float64,
L::Float64
)
if abs(x)<L/2 if abs(x)<L/2
return cos(pi*x/L)^2 return cos(pi*x/L)^2
else else
return 0. return 0.
end end
end end
# Fourier transform (in 3d)
@everywhere function hann_fourier(
k::Float64,
L::Float64
)
return L^2*4*pi^3/k*(((k*L)^3-4*k*L*pi^2)*cos(k*L/2)-2*(3*(k*L)^2-4*pi^2)*sin(k*L/2))/((k*L)^3-4*k*L*pi^2)^2
end
# normalized Gaussian (in 3d)
@everywhere function gaussian(
k::Float64,
L::Float64
)
return exp(-k^2/(2*L))/sqrt(8*pi^3*L^3)
end

36
src/interpolation.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -14,7 +14,11 @@
# linear interpolation: given vectors x,y, compute a linear interpolation for y(x0) # linear interpolation: given vectors x,y, compute a linear interpolation for y(x0)
# assume x is ordered # assume x is ordered
@everywhere function linear_interpolation(x0,x,y) @everywhere function linear_interpolation(
x0::Float64,
x::Array{Float64,1},
y::Array{Float64,1}
)
# if x0 is beyond all x's, then return the corresponding boundary value. # if x0 is beyond all x's, then return the corresponding boundary value.
if x0>x[length(x)] if x0>x[length(x)]
return y[length(y)] return y[length(y)]
@ -28,7 +32,12 @@
# interpolate # interpolate
return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i]) return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i])
end end
@everywhere function bracket(x0,x,a,b) @everywhere function bracket(
x0::Float64,
x::Array{Float64,1},
a::Int64,
b::Int64
)
i=floor(Int64,(a+b)/2) i=floor(Int64,(a+b)/2)
if x0<x[i] if x0<x[i]
return bracket(x0,x,a,i) return bracket(x0,x,a,i)
@ -41,15 +50,18 @@ end
# polynomial interpolation of a family of points # polynomial interpolation of a family of points
@everywhere function poly_interpolation(x,y) @everywhere function poly_interpolation(
x::Array{Float64,1},
y::Array{Float64,1}
)
# init for recursion # init for recursion
rec=Array{Polynomial{Float64}}(undef,length(x)) rec=Array{Polynomial{Float64},1}(undef,length(x))
for i in 1:length(x) for i in 1:length(x)
rec[i]=Polynomial([1.]) rec[i]=Polynomial([1.])
end end
# compute \prod (x-x_i) # compute \prod (x-x_i)
poly_interpolation_rec(rec,x,1,length(x)) poly_interpolation_rec(rec,x,1.,length(x))
# sum terms together # sum terms together
out=0. out=0.
@ -60,7 +72,12 @@ end
return out return out
end end
# recursive helper function # recursive helper function
@everywhere function poly_interpolation_rec(out,x,a,b) @everywhere function poly_interpolation_rec(
out::Array{Float64,1},
x::Array{Float64,1},
a::Float64,
b::Float64
)
if a==b if a==b
return return
end end
@ -91,7 +108,10 @@ end
return return
end end
## the following does the same, but has complexity N^2, the function above has N*log(N) ## the following does the same, but has complexity N^2, the function above has N*log(N)
#@everywhere function poly_interpolation(x,y) #@everywhere function poly_interpolation(
# x::Array{Float64,1},
# y::Array{Float64,1}
#)
# out=Polynomial([0.]) # out=Polynomial([0.])
# for i in 1:length(x) # for i in 1:length(x)
# prod=Polynomial([1.]) # prod=Polynomial([1.])

168
src/main.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -22,6 +22,8 @@ include("chebyshev.jl")
include("integration.jl") include("integration.jl")
include("interpolation.jl") include("interpolation.jl")
include("tools.jl") include("tools.jl")
include("multithread.jl")
include("optimization.jl")
include("potentials.jl") include("potentials.jl")
include("print.jl") include("print.jl")
include("easyeq.jl") include("easyeq.jl")
@ -36,9 +38,6 @@ function main()
rho=1e-6 rho=1e-6
e=1e-4 e=1e-4
# incrementally initialize Newton algorithm
nlrho_init=1
# potential # potential
v=k->v_exp(k,1.) v=k->v_exp(k,1.)
a0=a0_exp(1.) a0=a0_exp(1.)
@ -49,19 +48,29 @@ function main()
v_param_e=1. v_param_e=1.
# plot range when plotting in rho # plot range when plotting in rho
minlrho=-6 # linear
maxlrho=2 minrho=1e-6
maxrho=1e2
nrho=0
# logarithmic
minlrho=-6.
maxlrho=2.
nlrho=100 nlrho=100
rhos=Array{Float64}(undef,0) # list
rhos=Array{Float64,1}(undef,0)
# plot range when plotting in e # plot range when plotting in e
minle=-6 minle=-6.
maxle=2 maxle=2.
nle=100 nle=100
es=Array{Float64}(undef,0) es=Array{Float64,1}(undef,0)
# plot range when plotting in x # plot range when plotting in x
xmin=0 xmin=0.
xmax=100 xmax=100.
nx=100 nx=100
# plot range when plotting in k
kmin=0.
kmax=10.
nk=100
# cutoffs # cutoffs
tolerance=1e-11 tolerance=1e-11
@ -77,8 +86,8 @@ function main()
J=10 J=10
# starting rho from which to incrementally initialize Newton algorithm # starting rho from which to incrementally initialize Newton algorithm
# default must be set after reading rho, if not set explicitly minlrho_init=-6.
minlrho_init=nothing nlrho_init=0
# Hann window for Fourier transforms # Hann window for Fourier transforms
windowL=1e3 windowL=1e3
@ -98,6 +107,18 @@ function main()
anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.) anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.)
anyeq_approx=anyeq_bigeq_approx anyeq_approx=anyeq_bigeq_approx
# numerical approximations of derivatives
dx=1e-7
dk=1e-7
# initial guess for 2pt_max
x0=1.
k0=1.
# maximum step in 2pt_max
maxstep=Inf
# tolerance for max search
tolerance_max=Inf
# read cli arguments # read cli arguments
(params,potential,method,savefile,command)=read_args(ARGS) (params,potential,method,savefile,command)=read_args(ARGS)
@ -109,8 +130,8 @@ function main()
print(stderr,"error: could not read parameter '",param,"'.\n") print(stderr,"error: could not read parameter '",param,"'.\n")
exit(-1) exit(-1)
end end
lhs=terms[1] lhs=string(terms[1])
rhs=terms[2] rhs=string(terms[2])
if lhs=="rho" if lhs=="rho"
rho=parse(Float64,rhs) rho=parse(Float64,rhs)
elseif lhs=="minlrho_init" elseif lhs=="minlrho_init"
@ -133,6 +154,12 @@ function main()
maxlrho=parse(Float64,rhs) maxlrho=parse(Float64,rhs)
elseif lhs=="nlrho" elseif lhs=="nlrho"
nlrho=parse(Int64,rhs) nlrho=parse(Int64,rhs)
elseif lhs=="minrho"
minrho=parse(Float64,rhs)
elseif lhs=="maxrho"
maxrho=parse(Float64,rhs)
elseif lhs=="nrho"
nrho=parse(Int64,rhs)
elseif lhs=="es" elseif lhs=="es"
es=parse_list(rhs) es=parse_list(rhs)
elseif lhs=="minle" elseif lhs=="minle"
@ -147,6 +174,12 @@ function main()
xmax=parse(Float64,rhs) xmax=parse(Float64,rhs)
elseif lhs=="nx" elseif lhs=="nx"
nx=parse(Int64,rhs) nx=parse(Int64,rhs)
elseif lhs=="kmin"
kmin=parse(Float64,rhs)
elseif lhs=="kmax"
kmax=parse(Float64,rhs)
elseif lhs=="nk"
nk=parse(Int64,rhs)
elseif lhs=="P" elseif lhs=="P"
P=parse(Int64,rhs) P=parse(Int64,rhs)
elseif lhs=="N" elseif lhs=="N"
@ -155,6 +188,18 @@ function main()
J=parse(Int64,rhs) J=parse(Int64,rhs)
elseif lhs=="window_L" elseif lhs=="window_L"
windowL=parse(Float64,rhs) windowL=parse(Float64,rhs)
elseif lhs=="dx"
dx=parse(Float64,rhs)
elseif lhs=="x0"
x0=parse(Float64,rhs)
elseif lhs=="dk"
dk=parse(Float64,rhs)
elseif lhs=="k0"
k0=parse(Float64,rhs)
elseif lhs=="maxstep"
maxstep=parse(Float64,rhs)
elseif lhs=="tolerance_max"
tolerance_max=parse(Float64,rhs)
elseif lhs=="aK" elseif lhs=="aK"
anyeq_approx.aK=parse(Float64,rhs) anyeq_approx.aK=parse(Float64,rhs)
elseif lhs=="bK" elseif lhs=="bK"
@ -242,37 +287,48 @@ function main()
## set parameters ## set parameters
# rhos # rhos
if length(rhos)==0 if length(rhos)==0
rhos=Array{Float64}(undef,nlrho) # linear only if nrho is specified
for j in 0:nlrho-1 if nrho>0
rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j)) rhos=Array{Float64,1}(undef,nrho)
for j in 0:nrho-1
rhos[j+1]=(nrho==1 ? minrho : minrho+(maxrho-minrho)/(nrho-1)*j)
end
else
rhos=Array{Float64,1}(undef,nlrho)
for j in 0:nlrho-1
rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j))
end
end end
end end
# es # es
if length(es)==0 if length(es)==0
es=Array{Float64}(undef,nle) es=Array{Float64,1}(undef,nle)
for j in 0:nle-1 for j in 0:nle-1
es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j)) es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j))
end end
end end
# default minlrho_init
if (minlrho_init==nothing)
minlrho_init=log10(rho)
end
# splines # splines
taus=Array{Float64}(undef,J+1) taus=Array{Float64,1}(undef,J+1)
for j in 0:J for j in 0:J
taus[j+1]=-1+2*j/J taus[j+1]=-1+2*j/J
end end
# tolerance_max
if tolerance_max==Inf
tolerance_max=tolerance
end
## run command ## run command
if method=="easyeq" if command=="scattering_length"
@printf("% .15e\n",a0)
elseif method=="easyeq"
if command=="energy" if command=="energy"
easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
# e(rho) # e(rho)
elseif command=="energy_rho" elseif command=="energy_rho"
easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx) easyeq_energy_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx)
# u(k) # u(k)
elseif command=="uk" elseif command=="uk"
easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
@ -285,7 +341,28 @@ function main()
elseif command=="condensate_fraction" elseif command=="condensate_fraction"
easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
elseif command=="condensate_fraction_rho" elseif command=="condensate_fraction_rho"
easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx) easyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx)
# 2pt correlation
elseif command=="2pt"
easyeq_2pt(xmin,xmax,nx,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
# max of 2pt correlation
elseif command=="2pt_max"
easyeq_2pt_max(dx,x0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
# max of 2pt correlation as a function of rho
elseif command=="2pt_max_rho"
easyeq_2pt_max_rho(rhos,dx,x0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
# fourier transform of 2pt correlation
elseif command=="2pt_fourier"
easyeq_2pt_fourier(kmin,kmax,nk,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
# max of fourier transform of 2pt correlation
elseif command=="2pt_fourier_max"
easyeq_2pt_fourier_max(dk,k0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
# max of 2pt correlation as a function of rho
elseif command=="2pt_fourier_max_rho"
easyeq_2pt_fourier_max_rho(rhos,dk,k0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx)
# momentum distribution
elseif command=="momentum_distribution"
easyeq_momentum_distribution(kmin,kmax,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx)
else else
print(stderr,"unrecognized command '",command,"'.\n") print(stderr,"unrecognized command '",command,"'.\n")
exit(-1) exit(-1)
@ -316,7 +393,7 @@ function main()
anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
# e(rho) # e(rho)
elseif command=="energy_rho" elseif command=="energy_rho"
anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) anyeq_energy_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
elseif command=="energy_rho_init_prevrho" elseif command=="energy_rho_init_prevrho"
anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
elseif command=="energy_rho_init_nextrho" elseif command=="energy_rho_init_nextrho"
@ -331,12 +408,29 @@ function main()
elseif command=="condensate_fraction" elseif command=="condensate_fraction"
anyeq_condensate_fraction(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) anyeq_condensate_fraction(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
elseif command=="condensate_fraction_rho" elseif command=="condensate_fraction_rho"
anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) anyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
# momentum distribution # momentum distribution
elseif command=="momentum_distribution" elseif command=="momentum_distribution"
anyeq_momentum_distribution(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) anyeq_momentum_distribution(kmin,kmax,rho,minlrho_init,nlrho_init,taus,P,N,J,windowL,a0,v,maxiter,tolerance,anyeq_approx,savefile)
elseif command=="2pt" elseif command=="2pt"
anyeq_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile) anyeq_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
elseif command=="2pt_max"
anyeq_2pt_correlation_max(rho,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
elseif command=="2pt_max_rho"
anyeq_2pt_correlation_max_rho(rhos,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
elseif command=="2pt_fourier"
anyeq_2pt_correlation_fourier(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,kmin,kmax,nk,anyeq_approx,savefile)
elseif command=="2pt_fourier_test"
anyeq_2pt_correlation_fourier_test(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmax,kmin,kmax,nk,anyeq_approx,savefile)
elseif command=="2pt_fourier_max"
anyeq_2pt_correlation_fourier_max(rho,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
elseif command=="2pt_fourier_max_rho"
anyeq_2pt_correlation_fourier_max_rho(rhos,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile)
elseif command=="uncondensed_2pt"
anyeq_uncondensed_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
# compressibility
elseif command=="compressibility_rho"
anyeq_compressibility_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
else else
print(stderr,"unrecognized command: '",command,"'\n") print(stderr,"unrecognized command: '",command,"'\n")
exit(-1) exit(-1)
@ -360,9 +454,11 @@ function main()
end end
# parse a comma separated list as an array of Float64 # parse a comma separated list as an array of Float64
function parse_list(str) function parse_list(
str::String
)
elems=split(str,",") elems=split(str,",")
out=Array{Float64}(undef,length(elems)) out=Array{Float64,1}(undef,length(elems))
for i in 1:length(elems) for i in 1:length(elems)
out[i]=parse(Float64,elems[i]) out[i]=parse(Float64,elems[i])
end end
@ -370,7 +466,9 @@ function parse_list(str)
end end
# read cli arguments # read cli arguments
function read_args(ARGS) function read_args(
ARGS
)
# flag # flag
flag="" flag=""

16
src/multithread.jl Normal file
View File

@ -0,0 +1,16 @@
# split up 1...n among workers
function spawn_workers(n::Int64)
# number of workers
nw=nworkers()
# split jobs among workers
work=Array{Array{Int64,1},1}(undef,nw)
# init empty arrays
for p in 1:nw
work[p]=Int64[]
end
for i in 1:n
append!(work[(i-1)%nw+1],[i])
end
return work
end

94
src/optimization.jl Normal file
View File

@ -0,0 +1,94 @@
# gradient descent: find local minimum of function of one variable from initial guess
# numerically estimate the derivative
@everywhere function gradient_descent(
f::Function,
x0::Float64,
delta::Float64, # shift is delta*df
dx::Float64, # finite difference for numerical derivative evaluation
maxiter::Int64 # interrupt and fail after maxiter steps
)
counter=0
# init
x=x0
while counter<maxiter
# value at x and around
val=f(x)
valm=f(x-dx)
valp=f(x+dx)
# quit if minimum
if(val<valm && val<valp)
return(x,val)
end
# derivative
df=(valp-val)/dx
# step
x=x-delta*df
counter+=1
end
# fail
return(Inf,Inf)
end
# Newton algorithm to compute extrema
# numerically estimate the derivatives
@everywhere function newton_maximum(
f::Function,
x0::Float64,
dx::Float64, # finite difference for numerical derivative evaluation
maxiter::Int64,
tolerance::Float64,
maxstep::Float64 # maximal size of step
)
counter=0
# init
x=x0
while counter<maxiter
# value at x and around
val=f(x)
valm=f(x-dx)
valp=f(x+dx)
# derivative
dfp=(valp-val)/dx
dfm=(val-valm)/dx
# second derivative
ddf=(dfp-dfm)/dx
#@printf(stderr,"% .15e % .15e % .15e % .15e\n",x,val,dfp,ddf)
if abs(dfp/ddf)<tolerance
# check it is a local maximum
if ddf<0
return(x,val)
else
return(Inf,Inf)
end
end
# step
step=dfp/abs(ddf)
# step too large
if abs(step)>maxstep
step=maxstep*sign(step)
end
x=x+step
# fail if off to infinity
if x==Inf || x==-Inf
return(x,val)
end
counter+=1
end
# fail
return(Inf,Inf)
end

58
src/potentials.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,10 +13,15 @@
## limitations under the License. ## limitations under the License.
# exponential potential in 3 dimensions # exponential potential in 3 dimensions
@everywhere function v_exp(k,a) @everywhere function v_exp(
k::Float64,
a::Float64
)
return 8*pi/(1+k^2)^2*a return 8*pi/(1+k^2)^2*a
end end
@everywhere function a0_exp(a) @everywhere function a0_exp(
a::Float64
)
if a>0. if a>0.
return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a)) return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a))
elseif a<0. elseif a<0.
@ -27,15 +32,24 @@ end
end end
# exp(-x)-a*exp(-b*x) in 3 dimensions # exp(-x)-a*exp(-b*x) in 3 dimensions
@everywhere function v_expcry(k,a,b) @everywhere function v_expcry(
k::Float64,
a::Float64,
b::Float64
)
return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2)) return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2))
end end
@everywhere function a0_expcry(a,b) @everywhere function a0_expcry(
a::Float64,
b::Float64
)
return 1.21751642717932720441274114683413710125487579284827462 #ish return 1.21751642717932720441274114683413710125487579284827462 #ish
end end
# x^2*exp(-|x|) in 3 dimensions # x^2*exp(-|x|) in 3 dimensions
@everywhere function v_npt(k) @everywhere function v_npt(
k::Float64
)
return 96*pi*(1-k^2)/(1+k^2)^4 return 96*pi*(1-k^2)/(1+k^2)^4
end end
@everywhere function a0_npt() @everywhere function a0_npt()
@ -43,7 +57,9 @@ end
end end
# 1/(1+x^4/4) potential in 3 dimensions # 1/(1+x^4/4) potential in 3 dimensions
@everywhere function v_alg(k) @everywhere function v_alg(
k::Float64
)
if(k==0) if(k==0)
return 4*pi^2 return 4*pi^2
else else
@ -53,32 +69,50 @@ end
a0_alg=1. #ish a0_alg=1. #ish
# (1+a x^4)/(1+x^2)^4 potential in 3 dimensions # (1+a x^4)/(1+x^2)^4 potential in 3 dimensions
@everywhere function v_algwell(k) @everywhere function v_algwell(
k::Float64
)
a=4 a=4
return pi^2/24*exp(-k)*(a*(k^2-9*k+15)+k^2+3*k+3) return pi^2/24*exp(-k)*(a*(k^2-9*k+15)+k^2+3*k+3)
end end
a0_algwell=1. #ish a0_algwell=1. #ish
# potential corresponding to the exact solution c/(1+b^2x^2)^2 # potential corresponding to the exact solution c/(1+b^2x^2)^2
@everywhere function v_exact(k,b,c,e) @everywhere function v_exact(
k::Float64,
b::Float64,
c::Float64,
e::Float64
)
if k!=0 if k!=0
return 48*pi^2*((18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))*exp(-sqrt(1-sqrt(c))*k/b)+(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))*exp(-sqrt(1+sqrt(c))*k/b)+(1-k/b)/2*exp(-k/b)-c*e/b^2*(3*(9-c)*k/b+8*c)/(8*(9-c)^2)*exp(-2*k/b))/k return 48*pi^2*((18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))*exp(-sqrt(1-sqrt(c))*k/b)+(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))*exp(-sqrt(1+sqrt(c))*k/b)+(1-k/b)/2*exp(-k/b)-c*e/b^2*(3*(9-c)*k/b+8*c)/(8*(9-c)^2)*exp(-2*k/b))/k
else else
return 48*pi^2*(-sqrt(1-sqrt(c))/b*(18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))-sqrt(1+sqrt(c))/b*(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))-1/b-c*e/b^2*(27-19*c)/(8*(9-c)^2)) return 48*pi^2*(-sqrt(1-sqrt(c))/b*(18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))-sqrt(1+sqrt(c))/b*(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))-1/b-c*e/b^2*(27-19*c)/(8*(9-c)^2))
end end
end end
@everywhere function a0_exact(b,c,e) @everywhere function a0_exact(
b::Float64,
c::Float64,
e::Float64
)
return 1. #ish return 1. #ish
end end
# tent potential (convolution of soft sphere with itself): a*pi/12*(2*|x|/b-2)^2*(2*|x|/b+4) for |x|<b # tent potential (convolution of soft sphere with itself): a*pi/12*(2*|x|/b-2)^2*(2*|x|/b+4) for |x|<b
@everywhere function v_tent(k,a,b) @everywhere function v_tent(
k::Float64,
a::Float64,
b::Float64
)
if k!=0 if k!=0
return (b/2)^3*a*(4*pi*(sin(k*b/2)-k*b/2*cos(k*b/2))/(k*b/2)^3)^2 return (b/2)^3*a*(4*pi*(sin(k*b/2)-k*b/2*cos(k*b/2))/(k*b/2)^3)^2
else else
return (b/2)^3*a*(4*pi/3)^2 return (b/2)^3*a*(4*pi/3)^2
end end
end end
@everywhere function a0_tent(a,b) @everywhere function a0_tent(
a::Float64,
b::Float64
)
return b #ish return b #ish
end end

2
src/print.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.

68
src/simpleq-Kv.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,12 +13,27 @@
## limitations under the License. ## limitations under the License.
# Compute Kv=(-\Delta+v+4e(1-\rho u*))^{-1}v # Compute Kv=(-\Delta+v+4e(1-\rho u*))^{-1}v
function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) function anyeq_Kv(
minlrho::Float64,
nlrho::Int64,
taus::Array{Float64,1},
P::Int64,
N::Int64,
J::Int64,
rho::Float64,
a0::Float64,
v::Function,
maxiter::Int64,
tolerance::Float64,
xmin::Float64,
xmax::Float64,
nx::Int64
)
# init vectors # init vectors
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
# compute initial guess from medeq # compute initial guess from medeq
rhos=Array{Float64}(undef,nlrho) rhos=Array{Float64,1}(undef,nlrho)
for j in 0:nlrho-1 for j in 0:nlrho-1
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
end end
@ -44,7 +59,17 @@ function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,
end end
# Compute the condensate fraction for simpleq using Kv # Compute the condensate fraction for simpleq using Kv
function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance) function simpleq_Kv_condensate_fraction(
rhos::Array{Float64,1},
taus::Array{Float64,1},
P::Int64,
N::Int64,
J::Int64,
a0::Float64,
v::Function,
maxiter::Int64,
tolerance::Float64
)
# init vectors # init vectors
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
@ -67,19 +92,35 @@ function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
end end
# Compute the two-point correlation function for simpleq using Kv # Compute the two-point correlation function for simpleq using Kv
function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) function simpleq_Kv_2pt(
minlrho::Float64,
nlrho::Int64,
taus::Array{Float64,1},
P::Int64,
N::Int64,
J::Int64,
rho::Float64,
a0::Float64,
v::Function,
maxiter::Int64,
tolerance::Float64,
xmin::Float64,
xmax::Float64,
nx::Int64
)
# init vectors # init vectors
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
# compute initial guess from medeq # compute initial guess from medeq
rhos=Array{Float64}(undef,nlrho) rhos=Array{Float64,1}(undef,nlrho)
for j in 0:nlrho-1 for j in 0:nlrho-1
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
end end
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
u0=u0s[nlrho] u0=u0s[nlrho]
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,nothing,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)) Abar=Array{Float64,5}(undef,0,0,0,0,0)
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.))
# Kv in Fourier space # Kv in Fourier space
Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
@ -103,7 +144,18 @@ function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin
end end
# Kv # Kv
function simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) function simpleq_Kv_Kvk(
u::Array{Float64,1},
V::Array{Float64,1},
E::Float64,
rho::Float64,
Upsilon::Array{Array{Float64,1},1},
k::Array{Float64,1},
taus::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
# (-Delta+v+4e(1-\rho u*)) in Fourier space # (-Delta+v+4e(1-\rho u*)) in Fourier space
M=Array{Float64,2}(undef,N*J,N*J) M=Array{Float64,2}(undef,N*J,N*J)
for zetapp in 0:J-1 for zetapp in 0:J-1

486
src/simpleq-hardcore.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,28 +13,26 @@
## limitations under the License. ## limitations under the License.
# compute energy as a function of rho # compute energy as a function of rho
function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance) function simpleq_hardcore_energy_rho(
## spawn workers rhos::Array{Float64,1},
# number of workers taus::Array{Float64,1},
nw=nworkers() P::Int64,
# split jobs among workers N::Int64,
work=Array{Array{Int64,1},1}(undef,nw) J::Int64,
# init empty arrays maxiter::Int64,
for p in 1:nw tolerance::Float64
work[p]=zeros(0) )
end # spawn workers
for j in 1:length(rhos) work=spawn_workers(length(rhos))
append!(work[j%nw+1],j)
end
# initialize vectors # initialize vectors
(weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
# initial guess # initial guess
u0s=Array{Array{Float64}}(undef,length(rhos)) u0s=Array{Array{Float64,1}}(undef,length(rhos))
e0s=Array{Float64}(undef,length(rhos)) e0s=Array{Float64,1}(undef,length(rhos))
for j in 1:length(rhos) for j in 1:length(rhos)
u0s[j]=Array{Float64}(undef,N*J) u0s[j]=Array{Float64,1}(undef,N*J)
for i in 1:N*J for i in 1:N*J
u0s[j][i]=1/(1+r[i]^2)^2 u0s[j][i]=1/(1+r[i]^2)^2
end end
@ -43,17 +41,17 @@ function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
# save result from each task # save result from each task
us=Array{Array{Float64}}(undef,length(rhos)) us=Array{Array{Float64,1}}(undef,length(rhos))
es=Array{Float64}(undef,length(rhos)) es=Array{Float64,1}(undef,length(rhos))
err=Array{Float64}(undef,length(rhos)) err=Array{Float64,1}(undef,length(rhos))
count=0 count=0
# for each worker # for each worker
@sync for p in 1:nw @sync for p in 1:length(work)
# for each task # for each task
@async for j in work[p] @async for j in work[p]
count=count+1 count=count+1
if count>=nw if count>=length(work)
progress(count,length(rhos),10000) progress(count,length(rhos),10000)
end end
# run the task # run the task
@ -67,12 +65,20 @@ function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
end end
# compute u(x) # compute u(x)
function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance) function simpleq_hardcore_ux(
rho::Float64,
taus::Array{Float64,1},
P::Int64,
N::Int64,
J::Int64,
maxiter::Int64,
tolerance::Float64
)
# initialize vectors # initialize vectors
(weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
# initial guess # initial guess
u0=Array{Float64}(undef,N*J) u0=Array{Float64,1}(undef,N*J)
for i in 1:N*J for i in 1:N*J
u0[i]=1/(1+r[i]^2)^2 u0[i]=1/(1+r[i]^2)^2
end end
@ -86,28 +92,26 @@ function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
end end
# compute condensate fraction as a function of rho # compute condensate fraction as a function of rho
function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance) function simpleq_hardcore_condensate_fraction_rho(
## spawn workers rhos::Array{Float64,1},
# number of workers taus::Array{Float64,1},
nw=nworkers() P::Int64,
# split jobs among workers N::Int64,
work=Array{Array{Int64,1},1}(undef,nw) J::Int64,
# init empty arrays maxiter::Int64,
for p in 1:nw tolerance::Float64
work[p]=zeros(0) )
end # spawn workers
for j in 1:length(rhos) work=spawn_workers(length(rhos))
append!(work[j%nw+1],j)
end
# initialize vectors # initialize vectors
(weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
# initial guess # initial guess
u0s=Array{Array{Float64}}(undef,length(rhos)) u0s=Array{Array{Float64,1}}(undef,length(rhos))
e0s=Array{Float64}(undef,length(rhos)) e0s=Array{Float64,1}(undef,length(rhos))
for j in 1:length(rhos) for j in 1:length(rhos)
u0s[j]=Array{Float64}(undef,N*J) u0s[j]=Array{Float64,1}(undef,N*J)
for i in 1:N*J for i in 1:N*J
u0s[j][i]=1/(1+r[i]^2)^2 u0s[j][i]=1/(1+r[i]^2)^2
end end
@ -116,17 +120,17 @@ function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolera
# save result from each task # save result from each task
us=Array{Array{Float64}}(undef,length(rhos)) us=Array{Array{Float64,1}}(undef,length(rhos))
es=Array{Float64}(undef,length(rhos)) es=Array{Float64,1}(undef,length(rhos))
err=Array{Float64}(undef,length(rhos)) err=Array{Float64,1}(undef,length(rhos))
count=0 count=0
# for each worker # for each worker
@sync for p in 1:nw @sync for p in 1:length(work)
# for each task # for each task
@async for j in work[p] @async for j in work[p]
count=count+1 count=count+1
if count>=nw if count>=length(work)
progress(count,length(rhos),10000) progress(count,length(rhos),10000)
end end
# run the task # run the task
@ -142,13 +146,18 @@ end
# initialize computation # initialize computation
@everywhere function simpleq_hardcore_init(taus,P,N,J) @everywhere function simpleq_hardcore_init(
taus::Array{Float64,1},
P::Int64,
N::Int64,
J::Int64
)
# Gauss-Legendre weights # Gauss-Legendre weights
weights=gausslegendre(N) weights=gausslegendre(N)
weights_gL=gausslaguerre(N) weights_gL=gausslaguerre(N)
# r # r
r=Array{Float64}(undef,J*N) r=Array{Float64,1}(undef,J*N)
for zeta in 0:J-1 for zeta in 0:J-1
for j in 1:N for j in 1:N
xj=weights[1][j] xj=weights[1][j]
@ -164,9 +173,23 @@ end
end end
# compute u using chebyshev expansions # compute u using chebyshev expansions
@everywhere function simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) @everywhere function simpleq_hardcore_hatu(
u0::Array{Float64,1},
e0::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
maxiter::Int64,
tolerance::Float64
)
# init # init
vec=Array{Float64}(undef,J*N+1) vec=Array{Float64,1}(undef,J*N+1)
for i in 1:J*N for i in 1:J*N
vec[i]=u0[i] vec[i]=u0[i]
end end
@ -194,8 +217,20 @@ end
end end
# Xi # Xi
@everywhere function simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) @everywhere function simpleq_hardcore_Xi(
out=Array{Float64}(undef,J*N+1) u::Array{Float64,1},
e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64
)
out=Array{Float64,1}(undef,J*N+1)
FU=chebyshev(u,taus,weights,P,N,J,4) FU=chebyshev(u,taus,weights,P,N,J,4)
#D's #D's
@ -215,7 +250,19 @@ end
return out return out
end end
# DXi # DXi
@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) @everywhere function simpleq_hardcore_DXi(
u::Array{Float64,1},
e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64
)
out=Array{Float64,2}(undef,J*N+1,J*N+1) out=Array{Float64,2}(undef,J*N+1,J*N+1)
FU=chebyshev(u,taus,weights,P,N,J,4) FU=chebyshev(u,taus,weights,P,N,J,4)
@ -232,15 +279,15 @@ end
for zetapp in 0:J-1 for zetapp in 0:J-1
for n in 1:N for n in 1:N
one=zeros(Int64,J*N) one=zeros(Float64,J*N)
one[zetapp*N+n]=1 one[zetapp*N+n]=1.
Fone=chebyshev(one,taus,weights,P,N,J,4) Fone=chebyshev(one,taus,weights,P,N,J,4)
for i in 1:J*N for i in 1:J*N
# du/du # du/du
out[i,zetapp*N+n]=dot(Fone,d1[i])+2*dot(FU,d2[i]*Fone)-(zetapp*N+n==i ? 1 : 0) out[i,zetapp*N+n]=dot(Fone,d1[i])+2*dot(FU,d2[i]*Fone)-(zetapp*N+n==i ? 1 : 0)
# du/de # du/de
out[i,J*N+1]=(dsed0[i]+dot(FU,dsed1[i])+dot(FU,dsed2[i]*FU))/(2*sqrt(abs(e)))*(e>=0 ? 1 : -1) out[i,J*N+1]=(dsed0[i]+dot(FU,dsed1[i])+dot(FU,dsed2[i]*FU))/(2*sqrt(abs(e)))*(e>=0. ? 1. : -1.)
end end
# de/du # de/du
out[J*N+1,zetapp*N+n]=2*pi*rho* out[J*N+1,zetapp*N+n]=2*pi*rho*
@ -253,14 +300,26 @@ end
#de/de #de/de
out[J*N+1,J*N+1]=-1+2*pi*rho* out[J*N+1,J*N+1]=-1+2*pi*rho*
(2-dsedgamma0(e,rho,weights)-dot(FU,dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/denom/ (2-dsedgamma0(e,rho,weights)-dot(FU,dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/denom/
(2*sqrt(abs(e)))*(e>=0 ? 1 : -1) (2*sqrt(abs(e)))*(e>=0. ? 1. : -1.)
return out return out
end end
# dXi/dmu # dXi/dmu
@everywhere function simpleq_hardcore_dXidmu(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) @everywhere function simpleq_hardcore_dXidmu(
out=Array{Float64}(undef,J*N+1) u::Array{Float64,1},
e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64
)
out=Array{Float64,1}(undef,J*N+1)
FU=chebyshev(u,taus,weights,P,N,J,4) FU=chebyshev(u,taus,weights,P,N,J,4)
#D's #D's
@ -281,17 +340,37 @@ end
end end
# B's # B's
@everywhere function B0(r) @everywhere function B0(
r::Float64
)
return pi/12*(r-1)^2*(r+5) return pi/12*(r-1)^2*(r+5)
end end
@everywhere function B1(r,zeta,n,taus,T,weights,nu) @everywhere function B1(
r::Float64,
zeta::Int64,
n::Int64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu::Int64
)
return (taus[zeta+1]>=(2-r)/r || taus[zeta+2]<=-r/(r+2) ? 0 : return (taus[zeta+1]>=(2-r)/r || taus[zeta+2]<=-r/(r+2) ? 0 :
8*pi/(r+1)*integrate_legendre(tau-> 8*pi/(r+1)*integrate_legendre(tau->
(1-(r-(1-tau)/(1+tau))^2)/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],-r/(r+2)) (1-(r-(1-tau)/(1+tau))^2)/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],-r/(r+2))
,min(taus[zeta+2],(2-r)/r),weights) ,min(taus[zeta+2],(2-r)/r),weights)
) )
end end
@everywhere function B2(r,zeta,n,zetap,m,taus,T,weights,nu) @everywhere function B2(
r::Float64,
zeta::Int64,
n::Int64,
zetap::Int64,
m::Int64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
nu::Int64
)
return 32*pi/(r+1)*integrate_legendre(tau-> return 32*pi/(r+1)*integrate_legendre(tau->
1/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))* 1/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))*
(taus[zetap+1]>=alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau) || taus[zetap+2]<=alpham(1+r,tau) ? 0 : (taus[zetap+1]>=alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau) || taus[zetap+2]<=alpham(1+r,tau) ? 0 :
@ -303,30 +382,49 @@ end
end end
# D's # D's
@everywhere function D0(e,rho,r,weights,N,J) @everywhere function D0(
out=Array{Float64}(undef,J*N) e::Float64,
rho::Float64,
r::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=Array{Float64,1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ out[i]=exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s->
(s+1)*B0(s)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 (s+1)*B0(s)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2
,0,min(1,r[i]),weights)+ ,0.,min(1.,r[i]),weights)+
(r[i]>=1 ? 0 : (r[i]>=1 ? 0 :
rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_legendre(s-> rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_legendre(s->
(s+r[i]+1)*B0(s+r[i])*exp(-2*sqrt(abs(e))*s) (s+r[i]+1)*B0(s+r[i])*exp(-2*sqrt(abs(e))*s)
,0,1-r[i],weights) ,0.,1. -r[i],weights)
) )
end end
return out return out
end end
@everywhere function D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) @everywhere function D1(
out=Array{Array{Float64}}(undef,J*N) e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,1},1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=Array{Float64}(undef,(P+1)*J) out[i]=Array{Float64,1}(undef,(P+1)*J)
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[i][zeta*(P+1)+n+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> out[i][zeta*(P+1)+n+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s->
(s+1)*B1(s,zeta,n,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 (s+1)*B1(s,zeta,n,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2
,0,r[i],weights)+ ,0.,r[i],weights)+
rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s-> rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s->
(s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu) (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)
,2*sqrt(abs(e)),weights_gL) ,2*sqrt(abs(e)),weights_gL)
@ -336,7 +434,19 @@ end
return out return out
end end
@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) @everywhere function D2(
e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,2}}(undef,J*N) out=Array{Array{Float64,2}}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
@ -346,7 +456,7 @@ end
for m in 0:P for m in 0:P
out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s->
(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2
,0,r[i],weights)+ ,0.,r[i],weights)+
rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s-> rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s->
(s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu) (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)
,2*sqrt(abs(e)),weights_gL) ,2*sqrt(abs(e)),weights_gL)
@ -359,28 +469,47 @@ end
end end
# dD/d sqrt(abs(e))'s # dD/d sqrt(abs(e))'s
@everywhere function dseD0(e,rho,r,weights,N,J) @everywhere function dseD0(
out=Array{Float64}(undef,J*N) e::Float64,
rho::Float64,
r::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=Array{Float64,1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
rho/(r[i]+1)*integrate_legendre(s-> rho/(r[i]+1)*integrate_legendre(s->
(s+1)*B0(s)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) (s+1)*B0(s)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
,0,min(1,r[i]),weights)+ ,0.,min(1.,r[i]),weights)+
(r[i]>=1 ? 0 : (r[i]>=1 ? 0 :
rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights) rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0.,1. -r[i],weights)
) )
end end
return out return out
end end
@everywhere function dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) @everywhere function dseD1(
out=Array{Array{Float64}}(undef,J*N) e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,1},1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=Array{Float64}(undef,(P+1)*J) out[i]=Array{Float64,1}(undef,(P+1)*J)
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s-> out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s->
(s+1)*B1(s,zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
,0,r[i],weights)+ ,0.,r[i],weights)+
rho/(2*(r[i]+1))*integrate_laguerre(s-> rho/(2*(r[i]+1))*integrate_laguerre(s->
(s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
,2*sqrt(abs(e)),weights_gL) ,2*sqrt(abs(e)),weights_gL)
@ -389,7 +518,19 @@ end
end end
return out return out
end end
@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) @everywhere function dseD2(
e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,2}}(undef,J*N) out=Array{Array{Float64,2}}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
@ -399,7 +540,7 @@ end
for m in 0:P for m in 0:P
out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s-> out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s->
(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
,0,r[i],weights)+ ,0.,r[i],weights)+
rho/(2*(r[i]+1))*integrate_laguerre(s-> rho/(2*(r[i]+1))*integrate_laguerre(s->
(s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
,2*sqrt(abs(e)),weights_gL) ,2*sqrt(abs(e)),weights_gL)
@ -412,28 +553,47 @@ end
end end
# dD/d sqrt(abs(e+mu/2))'s # dD/d sqrt(abs(e+mu/2))'s
@everywhere function dsmuD0(e,rho,r,weights,N,J) @everywhere function dsmuD0(
out=Array{Float64}(undef,J*N) e::Float64,
rho::Float64,
r::Array{Float64,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
N::Int64,
J::Int64
)
out=Array{Float64,1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
rho/(r[i]+1)*integrate_legendre(s-> rho/(r[i]+1)*integrate_legendre(s->
(s+1)*B0(s)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) (s+1)*B0(s)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
,0,min(1,r[i]),weights)+ ,0.,min(1.,r[i]),weights)+
(r[i]>=1 ? 0 : (r[i]>=1 ? 0 :
rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights) rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0.,1. -r[i],weights)
) )
end end
return out return out
end end
@everywhere function dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) @everywhere function dsmuD1(
out=Array{Array{Float64}}(undef,J*N) e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,1},1}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=Array{Float64}(undef,(P+1)*J) out[i]=Array{Float64,1}(undef,(P+1)*J)
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s-> out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s->
(s+1)*B1(s,zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
,0,r[i],weights)+ ,0.,r[i],weights)+
rho/(2*(r[i]+1))*integrate_laguerre(s-> rho/(2*(r[i]+1))*integrate_laguerre(s->
(s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
,2*sqrt(abs(e)),weights_gL) ,2*sqrt(abs(e)),weights_gL)
@ -442,7 +602,19 @@ end
end end
return out return out
end end
@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) @everywhere function dsmuD2(
e::Float64,
rho::Float64,
r::Array{Float64,1},
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
N::Int64,
J::Int64,
nu::Int64
)
out=Array{Array{Float64,2}}(undef,J*N) out=Array{Array{Float64,2}}(undef,J*N)
for i in 1:J*N for i in 1:J*N
out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
@ -452,7 +624,7 @@ end
for m in 0:P for m in 0:P
out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s-> out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s->
(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2)
,0,r[i],weights)+ ,0.,r[i],weights)+
rho/(2*(r[i]+1))*integrate_laguerre(s-> rho/(2*(r[i]+1))*integrate_laguerre(s->
(s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i]))
,2*sqrt(abs(e)),weights_gL) ,2*sqrt(abs(e)),weights_gL)
@ -465,11 +637,25 @@ end
end end
# gamma's # gamma's
@everywhere function gamma0(e,rho,weights) @everywhere function gamma0(
return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights) e::Float64,
rho::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0.,1.,weights)
end end
@everywhere function gamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) @everywhere function gamma1(
out=Array{Float64}(undef,J*(P+1)) e::Float64,
rho::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,1}(undef,J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[zeta*(P+1)+n+1]=2*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) out[zeta*(P+1)+n+1]=2*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
@ -477,7 +663,17 @@ end
end end
return out return out
end end
@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) @everywhere function gamma2(
e::Float64,
rho::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
@ -492,11 +688,25 @@ end
end end
# dgamma/d sqrt(abs(e))'s # dgamma/d sqrt(abs(e))'s
@everywhere function dsedgamma0(e,rho,weights) @everywhere function dsedgamma0(
return 4*rho*sqrt(abs(e))*integrate_legendre(s->(s+1)*B0(s)*(1-sqrt(abs(e))*s)*exp(-2*sqrt(abs(e))*s),0,1,weights) e::Float64,
rho::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
return 4*rho*sqrt(abs(e))*integrate_legendre(s->(s+1)*B0(s)*(1-sqrt(abs(e))*s)*exp(-2*sqrt(abs(e))*s),0.,1.,weights)
end end
@everywhere function dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) @everywhere function dsedgamma1(
out=Array{Float64}(undef,J*(P+1)) e::Float64,
rho::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,1}(undef,J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[zeta*(P+1)+n+1]=4*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL) out[zeta*(P+1)+n+1]=4*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL)
@ -504,7 +714,17 @@ end
end end
return out return out
end end
@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) @everywhere function dsedgamma2(
e::Float64,
rho::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
@ -519,11 +739,25 @@ end
end end
# dgamma/d sqrt(e+mu/2)'s # dgamma/d sqrt(e+mu/2)'s
@everywhere function dsmudgamma0(e,rho,weights) @everywhere function dsmudgamma0(
return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights) e::Float64,
rho::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0.,1.,weights)
end end
@everywhere function dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) @everywhere function dsmudgamma1(
out=Array{Float64}(undef,J*(P+1)) e::Float64,
rho::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,1}(undef,J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[zeta*(P+1)+n+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) out[zeta*(P+1)+n+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
@ -531,7 +765,17 @@ end
end end
return out return out
end end
@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) @everywhere function dsmudgamma2(
e::Float64,
rho::Float64,
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
weights_gL::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
@ -546,25 +790,41 @@ end
end end
# \bar gamma's # \bar gamma's
@everywhere function gammabar0(weights) @everywhere function gammabar0(
return 4*pi*integrate_legendre(s->s^2*B0(s-1),0,1,weights) weights::Tuple{Array{Float64,1},Array{Float64,1}}
)
return 4*pi*integrate_legendre(s->s^2*B0(s-1),0.,1.,weights)
end end
@everywhere function gammabar1(taus,T,weights,P,J,nu) @everywhere function gammabar1(
out=Array{Float64}(undef,J*(P+1)) taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,1}(undef,J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
out[zeta*(P+1)+n+1]=4*pi*integrate_legendre(s->s^2*B1(s-1,zeta,n,taus,T,weights,nu),0,1,weights) out[zeta*(P+1)+n+1]=4*pi*integrate_legendre(s->s^2*B1(s-1,zeta,n,taus,T,weights,nu),0.,1.,weights)
end end
end end
return out return out
end end
@everywhere function gammabar2(taus,T,weights,P,J,nu) @everywhere function gammabar2(
taus::Array{Float64,1},
T::Array{Polynomial,1},
weights::Tuple{Array{Float64,1},Array{Float64,1}},
P::Int64,
J::Int64,
nu::Int64
)
out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
for zeta in 0:J-1 for zeta in 0:J-1
for n in 0:P for n in 0:P
for zetap in 0:J-1 for zetap in 0:J-1
for m in 0:P for m in 0:P
out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*pi*integrate_legendre(s->s^2*B2(s-1,zeta,n,zetap,m,taus,T,weights,nu),0,1,weights) out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*pi*integrate_legendre(s->s^2*B2(s-1,zeta,n,zetap,m,taus,T,weights,nu),0.,1.,weights)
end end
end end
end end

49
src/simpleq-iteration.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,7 +13,12 @@
## limitations under the License. ## limitations under the License.
# compute rho(e) using the iteration # compute rho(e) using the iteration
function simpleq_iteration_rho_e(es,order,v,maxiter) function simpleq_iteration_rho_e(
es::Array{Float64,1},
order::Int64,
v::Function,
maxiter::Int64
)
for j in 1:length(es) for j in 1:length(es)
(u,rho)=simpleq_iteration_hatun(es[j],order,v,maxiter) (u,rho)=simpleq_iteration_hatun(es[j],order,v,maxiter)
@printf("% .15e % .15e\n",es[j],real(rho[maxiter+1])) @printf("% .15e % .15e\n",es[j],real(rho[maxiter+1]))
@ -21,7 +26,15 @@ function simpleq_iteration_rho_e(es,order,v,maxiter)
end end
# compute u(x) using the iteration and print at every step # compute u(x) using the iteration and print at every step
function simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx) function simpleq_iteration_ux(
order::Int64,
e::Float64,
v::Function,
maxiter::Int64,
xmin::Float64,
xmax::Float64,
nx::Int64
)
(u,rho)=simpleq_iteration_hatun(e,order,v,maxiter) (u,rho)=simpleq_iteration_hatun(e,order,v,maxiter)
weights=gausslegendre(order) weights=gausslegendre(order)
@ -37,7 +50,12 @@ end
# \hat u_n # \hat u_n
function simpleq_iteration_hatun(e,order,v,maxiter) function simpleq_iteration_hatun(
e::Float64,
order::Int64,
v::Function,
maxiter::Int64
)
# gauss legendre weights # gauss legendre weights
weights=gausslegendre(order) weights=gausslegendre(order)
@ -46,7 +64,7 @@ function simpleq_iteration_hatun(e,order,v,maxiter)
(Eta,Eta0)=easyeq_init_H(weights,v) (Eta,Eta0)=easyeq_init_H(weights,v)
# init u and rho # init u and rho
u=Array{Array{Float64}}(undef,maxiter+1) u=Array{Array{Float64,1},1}(undef,maxiter+1)
u[1]=zeros(Float64,order) u[1]=zeros(Float64,order)
rho=zeros(Float64,maxiter+1) rho=zeros(Float64,maxiter+1)
@ -60,7 +78,11 @@ function simpleq_iteration_hatun(e,order,v,maxiter)
end end
# A # A
function simpleq_iteration_A(e,weights,Eta) function simpleq_iteration_A(
e::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}},
Eta::Array{Float64,1}
)
N=length(weights[1]) N=length(weights[1])
out=zeros(Float64,N,N) out=zeros(Float64,N,N)
for i in 1:N for i in 1:N
@ -75,7 +97,12 @@ function simpleq_iteration_A(e,weights,Eta)
end end
# b # b
function simpleq_iteration_b(u,e,rho,V) function simpleq_iteration_b(
u::Array{Float64,1},
e::Float64,
rho::Float64,
V::Array{Float64,1}
)
out=zeros(Float64,length(V)) out=zeros(Float64,length(V))
for i in 1:length(V) for i in 1:length(V)
out[i]=V[i]+2*e*rho*u[i]^2 out[i]=V[i]+2*e*rho*u[i]^2
@ -84,7 +111,13 @@ function simpleq_iteration_b(u,e,rho,V)
end end
# rho_n # rho_n
function simpleq_iteration_rhon(u,e,weights,V0,Eta0) function simpleq_iteration_rhon(
u::Array{Float64,1},
e::Float64,
weights::Tuple{Array{Float64,1},Array{Float64,1}},
V0::Float64,
Eta0::Float64
)
S=V0 S=V0
for i in 1:length(weights[1]) for i in 1:length(weights[1])
y=(weights[1][i]+1)/2 y=(weights[1][i]+1)/2

22
src/tools.jl
View File

@ -1,4 +1,4 @@
## Copyright 2021 Ian Jauslin ## Copyright 2021-2023 Ian Jauslin
## ##
## Licensed under the Apache License, Version 2.0 (the "License"); ## Licensed under the Apache License, Version 2.0 (the "License");
## you may not use this file except in compliance with the License. ## you may not use this file except in compliance with the License.
@ -13,7 +13,9 @@
## limitations under the License. ## limitations under the License.
# \Phi(x)=2*(1-sqrt(1-x))/x # \Phi(x)=2*(1-sqrt(1-x))/x
@everywhere function Phi(x) @everywhere function Phi(
x::Float64
)
if abs(x)>1e-5 if abs(x)>1e-5
return 2*(1-sqrt(abs(1-x)))/x return 2*(1-sqrt(abs(1-x)))/x
else else
@ -21,7 +23,9 @@
end end
end end
# \partial\Phi # \partial\Phi
@everywhere function dPhi(x) @everywhere function dPhi(
x::Float64
)
#if abs(x-1)<1e-5 #if abs(x-1)<1e-5
# @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x)) # @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x))
#end #end
@ -33,17 +37,25 @@ end
end end
# apply Phi to every element of a vector # apply Phi to every element of a vector
@everywhere function dotPhi(v) @everywhere function dotPhi(
v::Array{Float64,1}
)
out=zeros(Float64,length(v)) out=zeros(Float64,length(v))
for i in 1:length(v) for i in 1:length(v)
out[i]=Phi(v[i]) out[i]=Phi(v[i])
end end
return out return out
end end
@everywhere function dotdPhi(v) @everywhere function dotdPhi(
v::Array{Float64,1}
)
out=zeros(Float64,length(v)) out=zeros(Float64,length(v))
for i in 1:length(v) for i in 1:length(v)
out[i]=dPhi(v[i]) out[i]=dPhi(v[i])
end end
return out return out
end end
@everywhere function sinc(x::Float64)
return (x == 0 ? 1 : sin(x)/x)
end