Initial commit
This commit is contained in:
commit
e72af82c3e
202
LICENSE
Normal file
202
LICENSE
Normal file
@ -0,0 +1,202 @@
|
||||
|
||||
Apache License
|
||||
Version 2.0, January 2004
|
||||
http://www.apache.org/licenses/
|
||||
|
||||
TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
|
||||
|
||||
1. Definitions.
|
||||
|
||||
"License" shall mean the terms and conditions for use, reproduction,
|
||||
and distribution as defined by Sections 1 through 9 of this document.
|
||||
|
||||
"Licensor" shall mean the copyright owner or entity authorized by
|
||||
the copyright owner that is granting the License.
|
||||
|
||||
"Legal Entity" shall mean the union of the acting entity and all
|
||||
other entities that control, are controlled by, or are under common
|
||||
control with that entity. For the purposes of this definition,
|
||||
"control" means (i) the power, direct or indirect, to cause the
|
||||
direction or management of such entity, whether by contract or
|
||||
otherwise, or (ii) ownership of fifty percent (50%) or more of the
|
||||
outstanding shares, or (iii) beneficial ownership of such entity.
|
||||
|
||||
"You" (or "Your") shall mean an individual or Legal Entity
|
||||
exercising permissions granted by this License.
|
||||
|
||||
"Source" form shall mean the preferred form for making modifications,
|
||||
including but not limited to software source code, documentation
|
||||
source, and configuration files.
|
||||
|
||||
"Object" form shall mean any form resulting from mechanical
|
||||
transformation or translation of a Source form, including but
|
||||
not limited to compiled object code, generated documentation,
|
||||
and conversions to other media types.
|
||||
|
||||
"Work" shall mean the work of authorship, whether in Source or
|
||||
Object form, made available under the License, as indicated by a
|
||||
copyright notice that is included in or attached to the work
|
||||
(an example is provided in the Appendix below).
|
||||
|
||||
"Derivative Works" shall mean any work, whether in Source or Object
|
||||
form, that is based on (or derived from) the Work and for which the
|
||||
editorial revisions, annotations, elaborations, or other modifications
|
||||
represent, as a whole, an original work of authorship. For the purposes
|
||||
of this License, Derivative Works shall not include works that remain
|
||||
separable from, or merely link (or bind by name) to the interfaces of,
|
||||
the Work and Derivative Works thereof.
|
||||
|
||||
"Contribution" shall mean any work of authorship, including
|
||||
the original version of the Work and any modifications or additions
|
||||
to that Work or Derivative Works thereof, that is intentionally
|
||||
submitted to Licensor for inclusion in the Work by the copyright owner
|
||||
or by an individual or Legal Entity authorized to submit on behalf of
|
||||
the copyright owner. For the purposes of this definition, "submitted"
|
||||
means any form of electronic, verbal, or written communication sent
|
||||
to the Licensor or its representatives, including but not limited to
|
||||
communication on electronic mailing lists, source code control systems,
|
||||
and issue tracking systems that are managed by, or on behalf of, the
|
||||
Licensor for the purpose of discussing and improving the Work, but
|
||||
excluding communication that is conspicuously marked or otherwise
|
||||
designated in writing by the copyright owner as "Not a Contribution."
|
||||
|
||||
"Contributor" shall mean Licensor and any individual or Legal Entity
|
||||
on behalf of whom a Contribution has been received by Licensor and
|
||||
subsequently incorporated within the Work.
|
||||
|
||||
2. Grant of Copyright License. Subject to the terms and conditions of
|
||||
this License, each Contributor hereby grants to You a perpetual,
|
||||
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
|
||||
copyright license to reproduce, prepare Derivative Works of,
|
||||
publicly display, publicly perform, sublicense, and distribute the
|
||||
Work and such Derivative Works in Source or Object form.
|
||||
|
||||
3. Grant of Patent License. Subject to the terms and conditions of
|
||||
this License, each Contributor hereby grants to You a perpetual,
|
||||
worldwide, non-exclusive, no-charge, royalty-free, irrevocable
|
||||
(except as stated in this section) patent license to make, have made,
|
||||
use, offer to sell, sell, import, and otherwise transfer the Work,
|
||||
where such license applies only to those patent claims licensable
|
||||
by such Contributor that are necessarily infringed by their
|
||||
Contribution(s) alone or by combination of their Contribution(s)
|
||||
with the Work to which such Contribution(s) was submitted. If You
|
||||
institute patent litigation against any entity (including a
|
||||
cross-claim or counterclaim in a lawsuit) alleging that the Work
|
||||
or a Contribution incorporated within the Work constitutes direct
|
||||
or contributory patent infringement, then any patent licenses
|
||||
granted to You under this License for that Work shall terminate
|
||||
as of the date such litigation is filed.
|
||||
|
||||
4. Redistribution. You may reproduce and distribute copies of the
|
||||
Work or Derivative Works thereof in any medium, with or without
|
||||
modifications, and in Source or Object form, provided that You
|
||||
meet the following conditions:
|
||||
|
||||
(a) You must give any other recipients of the Work or
|
||||
Derivative Works a copy of this License; and
|
||||
|
||||
(b) You must cause any modified files to carry prominent notices
|
||||
stating that You changed the files; and
|
||||
|
||||
(c) You must retain, in the Source form of any Derivative Works
|
||||
that You distribute, all copyright, patent, trademark, and
|
||||
attribution notices from the Source form of the Work,
|
||||
excluding those notices that do not pertain to any part of
|
||||
the Derivative Works; and
|
||||
|
||||
(d) If the Work includes a "NOTICE" text file as part of its
|
||||
distribution, then any Derivative Works that You distribute must
|
||||
include a readable copy of the attribution notices contained
|
||||
within such NOTICE file, excluding those notices that do not
|
||||
pertain to any part of the Derivative Works, in at least one
|
||||
of the following places: within a NOTICE text file distributed
|
||||
as part of the Derivative Works; within the Source form or
|
||||
documentation, if provided along with the Derivative Works; or,
|
||||
within a display generated by the Derivative Works, if and
|
||||
wherever such third-party notices normally appear. The contents
|
||||
of the NOTICE file are for informational purposes only and
|
||||
do not modify the License. You may add Your own attribution
|
||||
notices within Derivative Works that You distribute, alongside
|
||||
or as an addendum to the NOTICE text from the Work, provided
|
||||
that such additional attribution notices cannot be construed
|
||||
as modifying the License.
|
||||
|
||||
You may add Your own copyright statement to Your modifications and
|
||||
may provide additional or different license terms and conditions
|
||||
for use, reproduction, or distribution of Your modifications, or
|
||||
for any such Derivative Works as a whole, provided Your use,
|
||||
reproduction, and distribution of the Work otherwise complies with
|
||||
the conditions stated in this License.
|
||||
|
||||
5. Submission of Contributions. Unless You explicitly state otherwise,
|
||||
any Contribution intentionally submitted for inclusion in the Work
|
||||
by You to the Licensor shall be under the terms and conditions of
|
||||
this License, without any additional terms or conditions.
|
||||
Notwithstanding the above, nothing herein shall supersede or modify
|
||||
the terms of any separate license agreement you may have executed
|
||||
with Licensor regarding such Contributions.
|
||||
|
||||
6. Trademarks. This License does not grant permission to use the trade
|
||||
names, trademarks, service marks, or product names of the Licensor,
|
||||
except as required for reasonable and customary use in describing the
|
||||
origin of the Work and reproducing the content of the NOTICE file.
|
||||
|
||||
7. Disclaimer of Warranty. Unless required by applicable law or
|
||||
agreed to in writing, Licensor provides the Work (and each
|
||||
Contributor provides its Contributions) on an "AS IS" BASIS,
|
||||
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
|
||||
implied, including, without limitation, any warranties or conditions
|
||||
of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A
|
||||
PARTICULAR PURPOSE. You are solely responsible for determining the
|
||||
appropriateness of using or redistributing the Work and assume any
|
||||
risks associated with Your exercise of permissions under this License.
|
||||
|
||||
8. Limitation of Liability. In no event and under no legal theory,
|
||||
whether in tort (including negligence), contract, or otherwise,
|
||||
unless required by applicable law (such as deliberate and grossly
|
||||
negligent acts) or agreed to in writing, shall any Contributor be
|
||||
liable to You for damages, including any direct, indirect, special,
|
||||
incidental, or consequential damages of any character arising as a
|
||||
result of this License or out of the use or inability to use the
|
||||
Work (including but not limited to damages for loss of goodwill,
|
||||
work stoppage, computer failure or malfunction, or any and all
|
||||
other commercial damages or losses), even if such Contributor
|
||||
has been advised of the possibility of such damages.
|
||||
|
||||
9. Accepting Warranty or Additional Liability. While redistributing
|
||||
the Work or Derivative Works thereof, You may choose to offer,
|
||||
and charge a fee for, acceptance of support, warranty, indemnity,
|
||||
or other liability obligations and/or rights consistent with this
|
||||
License. However, in accepting such obligations, You may act only
|
||||
on Your own behalf and on Your sole responsibility, not on behalf
|
||||
of any other Contributor, and only if You agree to indemnify,
|
||||
defend, and hold each Contributor harmless for any liability
|
||||
incurred by, or claims asserted against, such Contributor by reason
|
||||
of your accepting any such warranty or additional liability.
|
||||
|
||||
END OF TERMS AND CONDITIONS
|
||||
|
||||
APPENDIX: How to apply the Apache License to your work.
|
||||
|
||||
To apply the Apache License to your work, attach the following
|
||||
boilerplate notice, with the fields enclosed by brackets "[]"
|
||||
replaced with your own identifying information. (Don't include
|
||||
the brackets!) The text should be enclosed in the appropriate
|
||||
comment syntax for the file format. We also recommend that a
|
||||
file or class name and description of purpose be included on the
|
||||
same "printed page" as the copyright notice for easier
|
||||
identification within third-party archives.
|
||||
|
||||
Copyright [yyyy] [name of copyright owner]
|
||||
|
||||
Licensed under the Apache License, Version 2.0 (the "License");
|
||||
you may not use this file except in compliance with the License.
|
||||
You may obtain a copy of the License at
|
||||
|
||||
http://www.apache.org/licenses/LICENSE-2.0
|
||||
|
||||
Unless required by applicable law or agreed to in writing, software
|
||||
distributed under the License is distributed on an "AS IS" BASIS,
|
||||
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
See the License for the specific language governing permissions and
|
||||
limitations under the License.
|
||||
4
NOTICE
Normal file
4
NOTICE
Normal file
@ -0,0 +1,4 @@
|
||||
simplesolv:
|
||||
Licensed under the Apache 2.0 License (see LICENSE for details)
|
||||
Copyright Ian Jauslin 2021
|
||||
|
||||
102
README
Normal file
102
README
Normal file
@ -0,0 +1,102 @@
|
||||
simplesolv is a program to compute the solution of the Simplified approach to
|
||||
the Bose gas.
|
||||
|
||||
It is written in julia, and requires the julia interpreter to run.
|
||||
|
||||
The code is available in the 'src' directory.
|
||||
|
||||
The LaTeX source code for the documentation is located in the 'doc' directory.
|
||||
A pdf file is available in the source tarball, as well as at
|
||||
http://ian.jauslin.org/software/simplesolv/simplesolv-doc.pdf
|
||||
The documentation basic usage examples, an extensive reference of all commands,
|
||||
and detailed descriptions of the numerical algorithms used.
|
||||
|
||||
Any questions, concerns, or bug reports should be addressed to Ian Jauslin at
|
||||
ian.jauslin@rutgers.edu
|
||||
|
||||
simplesolv was written by Ian Jauslin, and is released under the Apache 2.0
|
||||
license. It is based on work done in collaboration with
|
||||
Eric A. Carlen
|
||||
Markus Holzmann
|
||||
Elliott H. Lieb
|
||||
published in the following papers:
|
||||
* E.A. Carlen, I. Jauslin, E.H. Lieb - 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,
|
||||
https://doi.org/10.2140/paa.2020.2.659
|
||||
https://arxiv.org/abs/1912.04987
|
||||
http://ian.jauslin.org/publications/19cjl
|
||||
|
||||
* 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
|
||||
Condensate Fraction, to appear in the SIAM Journal on Mathematical Analysis
|
||||
https://arxiv.org/abs/2010.13882
|
||||
http://ian.jauslin.org/publications/20cjl
|
||||
|
||||
* E.A. Carlen, M. Holzmann, I. Jauslin, E.H. Lieb - Simplified approach to
|
||||
the repulsive Bose gas from low to high densities and its numerical
|
||||
accuracy, Physical Review A, volume 103, issue 5, number 053309, 2021
|
||||
https://doi.org/10.1103/PhysRevA.103.053309
|
||||
https://arxiv.org/abs/2011.10869
|
||||
http://ian.jauslin.org/publications/20chjl
|
||||
|
||||
|
||||
|
||||
Dependencies:
|
||||
|
||||
* julia interpreter
|
||||
Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors: https://github.com/JuliaLang/julia/contributors
|
||||
https://julialang.org/
|
||||
MIT and GNU GPLv2 Licenses
|
||||
* julia packages:
|
||||
FastGaussQuadrature
|
||||
Alex Townsend
|
||||
https://github.com/JuliaApproximation/FastGaussQuadrature.jl
|
||||
MIT license
|
||||
Polynomials
|
||||
Jameson Nash and others
|
||||
https://github.com/JuliaMath/Polynomials.jl
|
||||
MIT license
|
||||
SpecialFunctions
|
||||
Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and others
|
||||
https://github.com/JuliaMath/SpecialFunctions.jl/
|
||||
MIT License
|
||||
|
||||
to typeset the documentation:
|
||||
* pdftex
|
||||
pdfTeX team, TeX Live Team, Hàn Thế Thành
|
||||
http://tug.org/applications/pdftex/
|
||||
GNU GPL
|
||||
* latex
|
||||
LaTeX3 Project
|
||||
https://www.latex-project.org/
|
||||
LaTeX Project Public License 1.3c
|
||||
* LaTeX packages:
|
||||
color
|
||||
David Carlisle, LaTeX3 Project
|
||||
https://ctan.org/pkg/color
|
||||
LaTeX Project Public License 1.3c
|
||||
marginnote
|
||||
Markus Kohm
|
||||
https://komascript.de/marginnote
|
||||
LaTeX Project Public License 1.3c
|
||||
amsfonts
|
||||
American Mathematical Society
|
||||
http://www.ams.org/tex/amsfonts.html
|
||||
SIL Open Font Licence
|
||||
hyperref
|
||||
Sebastian Rahtz, Heiko Oberdiek, LaTeX3 Project
|
||||
https://github.com/latex3/hyperref
|
||||
LaTeX Project Public License 1.3
|
||||
array
|
||||
LaTeX project, Frank Mittelbach
|
||||
https://ctan.org/pkg/array
|
||||
LaTeX Project Public License 1.3
|
||||
doublestroke
|
||||
Olaf Kummer
|
||||
https://ctan.org/pkg/doublestroke
|
||||
Custom Free license
|
||||
* optionally: GNU Make
|
||||
Richard Stallman, Roland McGrath, Paul Smith
|
||||
https://www.gnu.org/software/make/
|
||||
GNU GPLv3
|
||||
54
doc/Makefile
Normal file
54
doc/Makefile
Normal file
@ -0,0 +1,54 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
PROJECTNAME=$(basename $(wildcard *.tex))
|
||||
LIBS=$(notdir $(wildcard libs/*))
|
||||
|
||||
PDFS=$(addsuffix .pdf, $(PROJECTNAME))
|
||||
SYNCTEXS=$(addsuffix .synctex.gz, $(PROJECTNAME))
|
||||
|
||||
all: $(PROJECTNAME)
|
||||
|
||||
$(PROJECTNAME): $(LIBS)
|
||||
pdflatex -file-line-error $@.tex
|
||||
pdflatex -file-line-error $@.tex
|
||||
pdflatex -synctex=1 $@.tex
|
||||
|
||||
$(PROJECTNAME).aux: $(LIBS)
|
||||
pdflatex -file-line-error -draftmode $(PROJECTNAME).tex
|
||||
|
||||
|
||||
$(SYNCTEXS): $(LIBS)
|
||||
pdflatex -synctex=1 $(patsubst %.synctex.gz, %.tex, $@)
|
||||
|
||||
|
||||
libs: $(LIBS)
|
||||
|
||||
$(LIBS):
|
||||
ln -fs libs/$@ ./
|
||||
|
||||
|
||||
clean-aux:
|
||||
rm -f $(addsuffix .aux, $(PROJECTNAME))
|
||||
rm -f $(addsuffix .log, $(PROJECTNAME))
|
||||
rm -f $(addsuffix .out, $(PROJECTNAME))
|
||||
rm -f $(addsuffix .toc, $(PROJECTNAME))
|
||||
|
||||
clean-libs:
|
||||
rm -f $(LIBS)
|
||||
|
||||
clean-tex:
|
||||
rm -f $(PDFS) $(SYNCTEXS)
|
||||
|
||||
clean: clean-aux clean-tex clean-libs
|
||||
13
doc/bibliography/bibliography.tex
Normal file
13
doc/bibliography/bibliography.tex
Normal file
@ -0,0 +1,13 @@
|
||||
\bibitem[CHe21]{CHe21}E.A. Carlen, M. Holzmann, I. Jauslin, E.H. Lieb - {\it Simplified approach to the repulsive Bose gas from low to high densities and its numerical accuracy}, Physical Review A, volume~\-103, issue~\-5, number~\-053309, 2021,\par\penalty10000
|
||||
doi:{\tt\color{blue}\href{http://dx.doi.org/10.1103/PhysRevA.103.053309}{10.1103/PhysRevA.103.053309}}, arxiv:{\tt\color{blue}\href{http://arxiv.org/abs/2011.10869}{2011.10869}}.\par\medskip
|
||||
|
||||
\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
|
||||
|
||||
\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
|
||||
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[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
|
||||
|
||||
43
doc/libs/code.sty
Normal file
43
doc/libs/code.sty
Normal file
@ -0,0 +1,43 @@
|
||||
%% Copyright 2021 Ian Jauslin
|
||||
%%
|
||||
%% Licensed under the Apache License, Version 2.0 (the "License");
|
||||
%% you may not use this file except in compliance with the License.
|
||||
%% You may obtain a copy of the License at
|
||||
%%
|
||||
%% http://www.apache.org/licenses/LICENSE-2.0
|
||||
%%
|
||||
%% Unless required by applicable law or agreed to in writing, software
|
||||
%% distributed under the License is distributed on an "AS IS" BASIS,
|
||||
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
%% See the License for the specific language governing permissions and
|
||||
%% limitations under the License.
|
||||
|
||||
%%
|
||||
%% Code package:
|
||||
%% commands to typeset code
|
||||
%%
|
||||
|
||||
%% TeX format
|
||||
\NeedsTeXFormat{LaTeX2e}[1995/12/01]
|
||||
|
||||
%% package name
|
||||
\ProvidesPackage{code}[2021/03/15]
|
||||
|
||||
%% code environment
|
||||
\def\code{
|
||||
\par
|
||||
\leftskip10pt
|
||||
\bigskip
|
||||
\tt
|
||||
}
|
||||
\def\endcode{
|
||||
\rm
|
||||
\par
|
||||
\bigskip
|
||||
\leftskip0pt
|
||||
}
|
||||
|
||||
|
||||
%% end
|
||||
\endinput
|
||||
|
||||
55
doc/libs/dlmf.sty
Normal file
55
doc/libs/dlmf.sty
Normal file
@ -0,0 +1,55 @@
|
||||
%% Copyright 2021 Ian Jauslin
|
||||
%%
|
||||
%% Licensed under the Apache License, Version 2.0 (the "License");
|
||||
%% you may not use this file except in compliance with the License.
|
||||
%% You may obtain a copy of the License at
|
||||
%%
|
||||
%% http://www.apache.org/licenses/LICENSE-2.0
|
||||
%%
|
||||
%% Unless required by applicable law or agreed to in writing, software
|
||||
%% distributed under the License is distributed on an "AS IS" BASIS,
|
||||
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
%% See the License for the specific language governing permissions and
|
||||
%% limitations under the License.
|
||||
|
||||
%%
|
||||
%% DLMF package:
|
||||
%% cite equations from DLMF
|
||||
%%
|
||||
|
||||
%% TeX format
|
||||
\NeedsTeXFormat{LaTeX2e}[1995/12/01]
|
||||
|
||||
%% package name
|
||||
\ProvidesPackage{dlmf}[2020/05/01]
|
||||
|
||||
%% dependencies
|
||||
\RequirePackage{color}
|
||||
\RequirePackage{hyperref}
|
||||
|
||||
% get the first two numbers in (a.b.c);
|
||||
\def\@sectionnr(#1.#2.#3){#1.#2}
|
||||
|
||||
% get the last number in (a.b.c);
|
||||
\def\@eqnr(#1.#2.#3){#3}
|
||||
|
||||
% remove parentheses around argument
|
||||
\def\@cleanparentheses(#1){#1}
|
||||
|
||||
%% cite DLMF equation
|
||||
\def\dlmfcite#1#2{\leavevmode%
|
||||
\let\@dlmfcite@separator\@empty%
|
||||
\cite[%
|
||||
% loop over ',' separated list
|
||||
\@for\@dlmfcite:=#1\do{%
|
||||
% put commas between entries
|
||||
\@dlmfcite@separator\def\@dlmfcite@separator{,\ }%
|
||||
({\color{blue}\href{https://dlmf.nist.gov/\expandafter\@sectionnr\@dlmfcite\#E\expandafter\@eqnr\@dlmfcite}{\expandafter\@cleanparentheses\@dlmfcite}})%
|
||||
}%
|
||||
]{DLMF#2}%
|
||||
}
|
||||
|
||||
%% end
|
||||
\endinput
|
||||
|
||||
|
||||
682
doc/libs/ian.cls
Normal file
682
doc/libs/ian.cls
Normal file
@ -0,0 +1,682 @@
|
||||
%%
|
||||
%% Ian's class file
|
||||
%%
|
||||
|
||||
%% TeX format
|
||||
\NeedsTeXFormat{LaTeX2e}[1995/12/01]
|
||||
|
||||
%% class name
|
||||
\ProvidesClass{ian}[2017/09/29]
|
||||
|
||||
%% boolean to signal that this class is being used
|
||||
\newif\ifianclass
|
||||
|
||||
%% options
|
||||
% no section numbering in equations
|
||||
\DeclareOption{section_in_eq}{\sectionsineqtrue}
|
||||
\DeclareOption{section_in_fig}{\sectionsinfigtrue}
|
||||
\DeclareOption{section_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
|
||||
\DeclareOption{section_in_all}{\sectionsineqtrue\sectionsinfigtrue\PassOptionsToPackage{section_in_theo}{iantheo}}
|
||||
\DeclareOption{subsection_in_eq}{\subsectionsineqtrue}
|
||||
\DeclareOption{subsection_in_fig}{\subsectionsinfigtrue}
|
||||
\DeclareOption{subsection_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
|
||||
\DeclareOption{subsection_in_all}{\subsectionsineqtrue\subsectionsinfigtrue\PassOptionsToPackage{subsection_in_theo}{iantheo}}
|
||||
\DeclareOption{no_section_in_eq}{\sectionsineqfalse}
|
||||
\DeclareOption{no_section_in_fig}{\sectionsinfigfalse}
|
||||
\DeclareOption{no_section_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
|
||||
\DeclareOption{no_section_in_all}{\sectionsineqfalse\sectionsinfigfalse\PassOptionsToPackage{no_section_in_theo}{iantheo}}
|
||||
\DeclareOption{no_subsection_in_eq}{\subsectionsineqfalse}
|
||||
\DeclareOption{no_subsection_in_fig}{\subsectionsinfigfalse}
|
||||
\DeclareOption{no_subsection_in_theo}{\PassOptionsToPackage{\CurrentOption}{iantheo}}
|
||||
\DeclareOption{no_subsection_in_all}{\subsectionsineqfalse\subsectionsinfigfalse\PassOptionsToPackage{no_subsection_in_theo}{iantheo}}
|
||||
% reset point
|
||||
\DeclareOption{point_reset_at_section}{\PassOptionsToPackage{reset_at_section}{point}}
|
||||
\DeclareOption{point_no_reset_at_section}{\PassOptionsToPackage{no_reset_at_section}{point}}
|
||||
\DeclareOption{point_reset_at_theo}{\PassOptionsToPackage{reset_at_theo}{point}}
|
||||
\DeclareOption{point_no_reset_at_theo}{\PassOptionsToPackage{no_reset_at_theo}{point}}
|
||||
|
||||
\def\ian@defaultoptions{
|
||||
\ExecuteOptions{section_in_all, no_subsection_in_all}
|
||||
\ProcessOptions
|
||||
|
||||
%% required packages
|
||||
\RequirePackage{iantheo}
|
||||
\RequirePackage{point}
|
||||
\RequirePackage{color}
|
||||
\RequirePackage{marginnote}
|
||||
\RequirePackage{amssymb}
|
||||
\PassOptionsToPackage{hidelinks}{hyperref}
|
||||
\RequirePackage{hyperref}
|
||||
}
|
||||
|
||||
%% paper dimensions
|
||||
\setlength\paperheight{297mm}
|
||||
\setlength\paperwidth{210mm}
|
||||
|
||||
%% fonts
|
||||
\input{size11.clo}
|
||||
\DeclareOldFontCommand{\rm}{\normalfont\rmfamily}{\mathrm}
|
||||
\DeclareOldFontCommand{\sf}{\normalfont\sffamily}{\mathsf}
|
||||
\DeclareOldFontCommand{\tt}{\normalfont\ttfamily}{\mathtt}
|
||||
\DeclareOldFontCommand{\bf}{\normalfont\bfseries}{\mathbf}
|
||||
\DeclareOldFontCommand{\it}{\normalfont\itshape}{\mathit}
|
||||
|
||||
%% text dimensions
|
||||
\hoffset=-50pt
|
||||
\voffset=-72pt
|
||||
\textwidth=460pt
|
||||
\textheight=704pt
|
||||
|
||||
|
||||
%% remove default indentation
|
||||
\parindent=0pt
|
||||
%% indent command
|
||||
\def\indent{\hskip20pt}
|
||||
|
||||
%% something is wrong with \thepage, redefine it
|
||||
\gdef\thepage{\the\c@page}
|
||||
|
||||
%% array lines (to use the array environment)
|
||||
\setlength\arraycolsep{5\p@}
|
||||
\setlength\arrayrulewidth{.4\p@}
|
||||
|
||||
|
||||
%% correct vertical alignment at the end of a document
|
||||
\AtEndDocument{
|
||||
\vfill
|
||||
\eject
|
||||
}
|
||||
|
||||
|
||||
%% hyperlinks
|
||||
% hyperlinkcounter
|
||||
\newcounter{lncount}
|
||||
% hyperref anchor
|
||||
\def\hrefanchor{%
|
||||
\stepcounter{lncount}%
|
||||
\hypertarget{ln.\thelncount}{}%
|
||||
}
|
||||
|
||||
%% define a command and write it to aux file
|
||||
\def\outdef#1#2{%
|
||||
% define command%
|
||||
\expandafter\xdef\csname #1\endcsname{#2}%
|
||||
% hyperlink number%
|
||||
\expandafter\xdef\csname #1@hl\endcsname{\thelncount}%
|
||||
% write command to aux%
|
||||
\immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1\endcsname{\csname #1\endcsname}}%
|
||||
\immediate\write\@auxout{\noexpand\expandafter\noexpand\gdef\noexpand\csname #1@hl\endcsname{\thelncount}}%
|
||||
}
|
||||
|
||||
%% can call commands even when they are not defined
|
||||
\def\safe#1{%
|
||||
\ifdefined#1%
|
||||
#1%
|
||||
\else%
|
||||
{\color{red}\bf?}%
|
||||
\fi%
|
||||
}
|
||||
|
||||
%% define a label for the latest tag
|
||||
%% label defines a command containing the string stored in \tag
|
||||
\def\deflabel{
|
||||
\def\label##1{\expandafter\outdef{label@##1}{\safe\tag}}
|
||||
|
||||
\def\ref##1{%
|
||||
% check whether the label is defined (hyperlink runs into errors if this check is omitted)
|
||||
\ifcsname label@##1@hl\endcsname%
|
||||
\hyperlink{ln.\csname label@##1@hl\endcsname}{{\color{blue}\safe\csname label@##1\endcsname}}%
|
||||
\else%
|
||||
\ifcsname label@##1\endcsname%
|
||||
{\color{blue}\csname ##1\endcsname}%
|
||||
\else%
|
||||
{\bf ??}%
|
||||
\fi%
|
||||
\fi%
|
||||
}
|
||||
% manually specify label for ref
|
||||
\def\refname##1##2{%
|
||||
% check whether the label is defined (hyperlink runs into errors if this check is omitted)
|
||||
\ifcsname label@##1@hl\endcsname%
|
||||
\hyperlink{ln.\csname label@##1@hl\endcsname}{{\color{blue}##2}}%
|
||||
\else%
|
||||
??##2%
|
||||
\fi%
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
%% make a custom link at any given location in the document
|
||||
\def\makelink#1#2{%
|
||||
\hrefanchor%
|
||||
\outdef{label@#1}{#2}%
|
||||
}
|
||||
|
||||
|
||||
%% section command
|
||||
% counter
|
||||
\newcounter{sectioncount}
|
||||
% space before section
|
||||
\newlength\secskip
|
||||
\setlength\secskip{40pt}
|
||||
% a prefix to put before the section number, e.g. A for appendices
|
||||
\def\sectionprefix{}
|
||||
% define some lengths
|
||||
\newlength\secnumwidth
|
||||
\newlength\sectitlewidth
|
||||
\def\section#1{
|
||||
% reset counters
|
||||
\stepcounter{sectioncount}
|
||||
\setcounter{subsectioncount}{0}
|
||||
\ifsectionsineq
|
||||
\setcounter{seqcount}0
|
||||
\fi
|
||||
\ifsectionsinfig
|
||||
\setcounter{figcount}0
|
||||
\fi
|
||||
|
||||
% space before section (if not first)
|
||||
\ifnum\thesectioncount>1
|
||||
\vskip\secskip
|
||||
\penalty-1000
|
||||
\fi
|
||||
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\sectionprefix\thesectioncount}
|
||||
|
||||
% get widths
|
||||
\def\@secnum{{\bf\Large\sectionprefix\thesectioncount.\hskip10pt}}
|
||||
\settowidth\secnumwidth{\@secnum}
|
||||
\setlength\sectitlewidth\textwidth
|
||||
\addtolength\sectitlewidth{-\secnumwidth}
|
||||
|
||||
% print name
|
||||
\parbox{\textwidth}{
|
||||
\@secnum
|
||||
\parbox[t]{\sectitlewidth}{\Large\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@sec.\thesectioncount\endcsname{\thelncount}}
|
||||
\write\tocoutput{\noexpand\tocsection{#1}{\thepage}}
|
||||
\fi
|
||||
|
||||
%space
|
||||
\par\penalty10000
|
||||
\bigskip\penalty10000
|
||||
}
|
||||
|
||||
%% subsection
|
||||
% counter
|
||||
\newcounter{subsectioncount}
|
||||
% space before subsection
|
||||
\newlength\subsecskip
|
||||
\setlength\subsecskip{30pt}
|
||||
\def\subsection#1{
|
||||
% counters
|
||||
\stepcounter{subsectioncount}
|
||||
\setcounter{subsubsectioncount}{0}
|
||||
\ifsubsectionsineq
|
||||
\setcounter{seqcount}0
|
||||
\fi
|
||||
\ifsubsectionsinfig
|
||||
\setcounter{figcount}0
|
||||
\fi
|
||||
|
||||
% space before subsection (if not first)
|
||||
\ifnum\thesubsectioncount>1
|
||||
\vskip\subsecskip
|
||||
\penalty-500
|
||||
\fi
|
||||
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount}
|
||||
|
||||
% get widths
|
||||
\def\@secnum{{\bf\large\hskip.5cm\sectionprefix\thesectioncount.\thesubsectioncount.\hskip5pt}}
|
||||
\settowidth\secnumwidth{\@secnum}
|
||||
\setlength\sectitlewidth\textwidth
|
||||
\addtolength\sectitlewidth{-\secnumwidth}
|
||||
% print name
|
||||
\parbox{\textwidth}{
|
||||
\@secnum
|
||||
\parbox[t]{\sectitlewidth}{\large\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@subsec.\thesectioncount.\thesubsectioncount\endcsname{\thelncount}}
|
||||
\write\tocoutput{\noexpand\tocsubsection{#1}{\thepage}}
|
||||
\fi
|
||||
|
||||
% space
|
||||
\par\penalty10000
|
||||
\medskip\penalty10000
|
||||
}
|
||||
|
||||
%% subsubsection
|
||||
% counter
|
||||
\newcounter{subsubsectioncount}
|
||||
% space before subsubsection
|
||||
\newlength\subsubsecskip
|
||||
\setlength\subsubsecskip{20pt}
|
||||
\def\subsubsection#1{
|
||||
% counters
|
||||
\stepcounter{subsubsectioncount}
|
||||
|
||||
% space before subsubsection (if not first)
|
||||
\ifnum\thesubsubsectioncount>1
|
||||
\vskip\subsubsecskip
|
||||
\penalty-500
|
||||
\fi
|
||||
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount}
|
||||
|
||||
% get widths
|
||||
\def\@secnum{{\bf\hskip1.cm\sectionprefix\thesectioncount.\thesubsectioncount.\thesubsubsectioncount.\hskip5pt}}
|
||||
\settowidth\secnumwidth{\@secnum}
|
||||
\setlength\sectitlewidth\textwidth
|
||||
\addtolength\sectitlewidth{-\secnumwidth}
|
||||
% print name
|
||||
\parbox{\textwidth}{
|
||||
\@secnum
|
||||
\parbox[t]{\sectitlewidth}{\large\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@subsubsec.\thesectioncount.\thesubsectioncount.\thesubsubsectioncount\endcsname{\thelncount}}
|
||||
\write\tocoutput{\noexpand\tocsubsubsection{#1}{\thepage}}
|
||||
\fi
|
||||
|
||||
% space
|
||||
\par\penalty10000
|
||||
\medskip\penalty10000
|
||||
}
|
||||
|
||||
%% itemize
|
||||
\newlength\itemizeskip
|
||||
% left margin for items
|
||||
\setlength\itemizeskip{20pt}
|
||||
\newlength\itemizeseparator
|
||||
% space between the item symbol and the text
|
||||
\setlength\itemizeseparator{5pt}
|
||||
% penalty preceding an itemize
|
||||
\newcount\itemizepenalty
|
||||
\itemizepenalty=0
|
||||
% counter counting the itemize level
|
||||
\newcounter{itemizecount}
|
||||
|
||||
% item symbol
|
||||
\def\itemizept#1{
|
||||
\ifnum#1=1
|
||||
\textbullet
|
||||
\else
|
||||
$\scriptstyle\blacktriangleright$
|
||||
\fi
|
||||
}
|
||||
|
||||
|
||||
\newlength\current@itemizeskip
|
||||
\setlength\current@itemizeskip{0pt}
|
||||
\def\itemize{%
|
||||
\par\expandafter\penalty\the\itemizepenalty\medskip\expandafter\penalty\the\itemizepenalty%
|
||||
\addtocounter{itemizecount}{1}%
|
||||
\addtolength\current@itemizeskip{\itemizeskip}%
|
||||
\leftskip\current@itemizeskip%
|
||||
}
|
||||
\def\enditemize{%
|
||||
\addtocounter{itemizecount}{-1}%
|
||||
\addtolength\current@itemizeskip{-\itemizeskip}%
|
||||
\par\expandafter\penalty\the\itemizepenalty\leftskip\current@itemizeskip%
|
||||
\medskip\expandafter\penalty\the\itemizepenalty%
|
||||
}
|
||||
|
||||
% item, with optional argument to specify the item point
|
||||
% @itemarg is set to true when there is an optional argument
|
||||
\newif\if@itemarg
|
||||
\def\item{%
|
||||
% check whether there is an optional argument (if there is none, add on empty '[]')
|
||||
\@ifnextchar [{\@itemargtrue\@itemx}{\@itemargfalse\@itemx[]}%
|
||||
}
|
||||
\newlength\itempt@total
|
||||
\def\@itemx[#1]{
|
||||
\if@itemarg
|
||||
\settowidth\itempt@total{#1}
|
||||
\else
|
||||
\settowidth\itempt@total{\itemizept\theitemizecount}
|
||||
\fi
|
||||
\addtolength\itempt@total{\itemizeseparator}
|
||||
\par
|
||||
\medskip
|
||||
\if@itemarg
|
||||
\hskip-\itempt@total#1\hskip\itemizeseparator
|
||||
\else
|
||||
\hskip-\itempt@total\itemizept\theitemizecount\hskip\itemizeseparator
|
||||
\fi
|
||||
}
|
||||
|
||||
%% prevent page breaks after itemize
|
||||
\newcount\previtemizepenalty
|
||||
\def\nopagebreakafteritemize{
|
||||
\previtemizepenalty=\itemizepenalty
|
||||
\itemizepenalty=10000
|
||||
}
|
||||
%% back to previous value
|
||||
\def\restorepagebreakafteritemize{
|
||||
\itemizepenalty=\previtemizepenalty
|
||||
}
|
||||
|
||||
%% enumerate
|
||||
\newcounter{enumerate@count}
|
||||
\def\enumerate{
|
||||
\setcounter{enumerate@count}0
|
||||
\let\olditem\item
|
||||
\let\olditemizept\itemizept
|
||||
\def\item{
|
||||
% counter
|
||||
\stepcounter{enumerate@count}
|
||||
% set header
|
||||
\def\itemizept{\theenumerate@count.}
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\theenumerate@count}
|
||||
\olditem
|
||||
}
|
||||
\itemize
|
||||
}
|
||||
\def\endenumerate{
|
||||
\enditemize
|
||||
\let\item\olditem
|
||||
\let\itemizept\olditemizept
|
||||
}
|
||||
|
||||
|
||||
%% equation numbering
|
||||
% counter
|
||||
\newcounter{seqcount}
|
||||
% booleans (write section or subsection in equation number)
|
||||
\newif\ifsectionsineq
|
||||
\newif\ifsubsectionsineq
|
||||
\def\seqcount{
|
||||
\stepcounter{seqcount}
|
||||
% the output
|
||||
\edef\seqformat{\theseqcount}
|
||||
% add subsection number
|
||||
\ifsubsectionsineq
|
||||
\let\tmp\seqformat
|
||||
\edef\seqformat{\thesubsectioncount.\tmp}
|
||||
\fi
|
||||
% add section number
|
||||
\ifsectionsineq
|
||||
\let\tmp\seqformat
|
||||
\edef\seqformat{\sectionprefix\thesectioncount.\tmp}
|
||||
\fi
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\seqformat}
|
||||
% write number
|
||||
\marginnote{\hfill(\seqformat)}
|
||||
}
|
||||
%% equation environment compatibility
|
||||
\def\equation{\hrefanchor$$\seqcount}
|
||||
\def\endequation{$$\@ignoretrue}
|
||||
|
||||
|
||||
%% figures
|
||||
% counter
|
||||
\newcounter{figcount}
|
||||
% booleans (write section or subsection in equation number)
|
||||
\newif\ifsectionsinfig
|
||||
\newif\ifsubsectionsinfig
|
||||
% width of figures
|
||||
\newlength\figwidth
|
||||
\setlength\figwidth\textwidth
|
||||
\addtolength\figwidth{-2.5cm}
|
||||
% caption
|
||||
\def\defcaption{
|
||||
\long\def\caption##1{
|
||||
\stepcounter{figcount}
|
||||
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
|
||||
% the number of the figure
|
||||
\edef\figformat{\thefigcount}
|
||||
% add subsection number
|
||||
\ifsubsectionsinfig
|
||||
\let\tmp\figformat
|
||||
\edef\figformat{\thesubsectioncount.\tmp}
|
||||
\fi
|
||||
% add section number
|
||||
\ifsectionsinfig
|
||||
\let\tmp\figformat
|
||||
\edef\figformat{\sectionprefix\thesectioncount.\tmp}
|
||||
\fi
|
||||
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\figformat}
|
||||
|
||||
% write
|
||||
\hfil fig \figformat: \parbox[t]{\figwidth}{\leavevmode\small##1}
|
||||
|
||||
% space
|
||||
\par\bigskip
|
||||
}
|
||||
}
|
||||
%% short caption: centered
|
||||
\def\captionshort#1{
|
||||
\stepcounter{figcount}
|
||||
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
|
||||
% the number of the figure
|
||||
\edef\figformat{\thefigcount}
|
||||
% add section number
|
||||
\ifsectionsinfig
|
||||
\let\tmp\figformat
|
||||
\edef\figformat{\sectionprefix\thesectioncount.\tmp}
|
||||
\fi
|
||||
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\figformat}
|
||||
|
||||
% write
|
||||
\hfil fig \figformat: {\small#1}
|
||||
|
||||
%space
|
||||
\par\bigskip
|
||||
}
|
||||
|
||||
%% environment
|
||||
\def\figure{
|
||||
\par
|
||||
\vfil\penalty100\vfilneg
|
||||
\bigskip
|
||||
}
|
||||
\def\endfigure{
|
||||
\par
|
||||
\bigskip
|
||||
}
|
||||
|
||||
|
||||
%% start appendices
|
||||
\def\appendix{
|
||||
\vfill
|
||||
\pagebreak
|
||||
|
||||
% counter
|
||||
\setcounter{sectioncount}0
|
||||
|
||||
% prefix
|
||||
\def\sectionprefix{A}
|
||||
|
||||
% write
|
||||
{\bf \LARGE Appendices}\par\penalty10000\bigskip\penalty10000
|
||||
|
||||
% add a mention in the table of contents
|
||||
\iftoc
|
||||
\immediate\write\tocoutput{\noexpand\tocappendices}\penalty10000
|
||||
\fi
|
||||
|
||||
%% uncomment for new page for each appendix
|
||||
%\def\seqskip{\vfill\pagebreak}
|
||||
}
|
||||
|
||||
|
||||
%% bibliography
|
||||
% size of header
|
||||
\newlength\bibheader
|
||||
\def\thebibliography#1{
|
||||
\hrefanchor
|
||||
|
||||
% add a mention in the table of contents
|
||||
\iftoc
|
||||
% save lncount in aux variable which is written to toc
|
||||
\immediate\write\tocoutput{\noexpand\expandafter\noexpand\edef\noexpand\csname toc@references\endcsname{\thelncount}}
|
||||
\write\tocoutput{\noexpand\tocreferences{\thepage}}\penalty10000
|
||||
\fi
|
||||
|
||||
% write
|
||||
{\bf \LARGE References}\par\penalty10000\bigskip\penalty10000
|
||||
% width of header
|
||||
\settowidth\bibheader{[#1]}
|
||||
\leftskip\bibheader
|
||||
}
|
||||
% end environment
|
||||
\def\endthebibliography{
|
||||
\par\leftskip0pt
|
||||
}
|
||||
|
||||
%% bibitem command
|
||||
\def\bibitem[#1]#2{%
|
||||
\hrefanchor%
|
||||
\outdef{label@cite#2}{#1}%
|
||||
\hskip-\bibheader%
|
||||
\makebox[\bibheader]{\cite{#2}\hfill}%
|
||||
}
|
||||
|
||||
%% cite command
|
||||
% @tempswa is set to true when there is an optional argument
|
||||
\newif\@tempswa
|
||||
\def\cite{%
|
||||
% check whether there is an optional argument (if there is none, add on empty '[]')
|
||||
\@ifnextchar [{\@tempswatrue\@citex}{\@tempswafalse\@citex[]}%
|
||||
}
|
||||
% command with optional argument
|
||||
\def\@citex[#1]#2{\leavevmode%
|
||||
% initialize loop
|
||||
\let\@cite@separator\@empty%
|
||||
% format
|
||||
\@cite{%
|
||||
% loop over ',' separated list
|
||||
\@for\@cite@:=#2\do{%
|
||||
% text to add at each iteration of the loop (separator between citations)
|
||||
\@cite@separator\def\@cite@separator{,\ }%
|
||||
% add entry to citelist
|
||||
\@writecitation{\@cite@}%
|
||||
\ref{cite\@cite@}%
|
||||
}%
|
||||
}%
|
||||
% add optional argument text (as an argument to '\@cite')
|
||||
{#1}%
|
||||
}
|
||||
\def\@cite#1#2{%
|
||||
[#1\if@tempswa , #2\fi]%
|
||||
}
|
||||
%% add entry to citelist after checking it has not already been added
|
||||
\def\@writecitation#1{%
|
||||
\ifcsname if#1cited\endcsname%
|
||||
\else%
|
||||
\expandafter\newif\csname if#1cited\endcsname%
|
||||
\immediate\write\@auxout{\string\citation{#1}}%
|
||||
\fi%
|
||||
}
|
||||
|
||||
%% table of contents
|
||||
% boolean
|
||||
\newif\iftoc
|
||||
\def\tableofcontents{
|
||||
{\bf \large Table of contents:}\par\penalty10000\bigskip\penalty10000
|
||||
|
||||
% copy content from file
|
||||
\IfFileExists{\jobname.toc}{\input{\jobname.toc}}{{\tt error: table of contents missing}}
|
||||
|
||||
% open new toc
|
||||
\newwrite\tocoutput
|
||||
\immediate\openout\tocoutput=\jobname.toc
|
||||
|
||||
\toctrue
|
||||
}
|
||||
%% close file
|
||||
\AtEndDocument{
|
||||
% close toc
|
||||
\iftoc
|
||||
\immediate\closeout\tocoutput
|
||||
\fi
|
||||
}
|
||||
|
||||
|
||||
%% fill line with dots
|
||||
\def\leaderfill{\leaders\hbox to 1em {\hss. \hss}\hfill}
|
||||
|
||||
%% same as sectionprefix
|
||||
\def\tocsectionprefix{}
|
||||
|
||||
%% toc formats
|
||||
\newcounter{tocsectioncount}
|
||||
\def\tocsection #1#2{
|
||||
\stepcounter{tocsectioncount}
|
||||
\setcounter{tocsubsectioncount}{0}
|
||||
\setcounter{tocsubsubsectioncount}{0}
|
||||
% write
|
||||
\smallskip\hyperlink{ln.\csname toc@sec.\thetocsectioncount\endcsname}{{\bf \tocsectionprefix\thetocsectioncount}.\hskip5pt {\color{blue}#1}\leaderfill#2}\par
|
||||
}
|
||||
\newcounter{tocsubsectioncount}
|
||||
\def\tocsubsection #1#2{
|
||||
\stepcounter{tocsubsectioncount}
|
||||
\setcounter{tocsubsubsectioncount}{0}
|
||||
% write
|
||||
{\hskip10pt\hyperlink{ln.\csname toc@subsec.\thetocsectioncount.\thetocsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
|
||||
}
|
||||
\newcounter{tocsubsubsectioncount}
|
||||
\def\tocsubsubsection #1#2{
|
||||
\stepcounter{tocsubsubsectioncount}
|
||||
% write
|
||||
{\hskip20pt\hyperlink{ln.\csname toc@subsubsec.\thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount\endcsname}{{\bf \thetocsectioncount.\thetocsubsectioncount.\thetocsubsubsectioncount}.\hskip5pt {\color{blue}\small #1}\leaderfill#2}}\par
|
||||
}
|
||||
\def\tocappendices{
|
||||
\medskip
|
||||
\setcounter{tocsectioncount}0
|
||||
{\bf Appendices}\par
|
||||
\smallskip
|
||||
\def\tocsectionprefix{A}
|
||||
}
|
||||
\def\tocreferences#1{
|
||||
\medskip
|
||||
{\hyperlink{ln.\csname toc@references\endcsname}{{\color{blue}\bf References}\leaderfill#1}}\par
|
||||
\smallskip
|
||||
}
|
||||
|
||||
|
||||
%% definitions that must be loaded at begin document
|
||||
\let\ian@olddocument\document
|
||||
\def\document{
|
||||
\ian@olddocument
|
||||
|
||||
\deflabel
|
||||
\defcaption
|
||||
}
|
||||
|
||||
%% end
|
||||
\ian@defaultoptions
|
||||
\endinput
|
||||
176
doc/libs/iantheo.sty
Normal file
176
doc/libs/iantheo.sty
Normal file
@ -0,0 +1,176 @@
|
||||
%% Copyright 2021 Ian Jauslin
|
||||
%%
|
||||
%% Licensed under the Apache License, Version 2.0 (the "License");
|
||||
%% you may not use this file except in compliance with the License.
|
||||
%% You may obtain a copy of the License at
|
||||
%%
|
||||
%% http://www.apache.org/licenses/LICENSE-2.0
|
||||
%%
|
||||
%% Unless required by applicable law or agreed to in writing, software
|
||||
%% distributed under the License is distributed on an "AS IS" BASIS,
|
||||
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
%% See the License for the specific language governing permissions and
|
||||
%% limitations under the License.
|
||||
|
||||
%%
|
||||
%% iantheorem package:
|
||||
%% Ian's customized theorem command
|
||||
%%
|
||||
|
||||
%% boolean to signal that this package was loaded
|
||||
\newif\ifiantheo
|
||||
|
||||
%% TeX format
|
||||
\NeedsTeXFormat{LaTeX2e}[1995/12/01]
|
||||
|
||||
%% package name
|
||||
\ProvidesPackage{iantheo}[2016/11/10]
|
||||
|
||||
%% options
|
||||
\newif\ifsectionintheo
|
||||
\DeclareOption{section_in_theo}{\sectionintheotrue}
|
||||
\DeclareOption{no_section_in_theo}{\sectionintheofalse}
|
||||
\newif\ifsubsectionintheo
|
||||
\DeclareOption{subsection_in_theo}{\subsectionintheotrue}
|
||||
\DeclareOption{no_subsection_in_theo}{\subsectionintheofalse}
|
||||
|
||||
\def\iantheo@defaultoptions{
|
||||
\ExecuteOptions{section_in_theo, no_subsection_in_theo}
|
||||
\ProcessOptions
|
||||
|
||||
%%% reset at every new section
|
||||
\ifsectionintheo
|
||||
\let\iantheo@oldsection\section
|
||||
\gdef\section{\setcounter{theocount}{0}\iantheo@oldsection}
|
||||
\fi
|
||||
|
||||
%% reset at every new subsection
|
||||
\ifsubsectionintheo
|
||||
\let\iantheo@oldsubsection\subsection
|
||||
\gdef\subsection{\setcounter{theocount}{0}\iantheo@oldsubsection}
|
||||
\fi
|
||||
}
|
||||
|
||||
|
||||
%% delimiters
|
||||
\def\delimtitle#1{
|
||||
\par%
|
||||
\leavevmode%
|
||||
\raise.3em\hbox to\hsize{%
|
||||
\lower0.3em\hbox{\vrule height0.3em}%
|
||||
\hrulefill%
|
||||
\ \lower.3em\hbox{#1}\ %
|
||||
\hrulefill%
|
||||
\lower0.3em\hbox{\vrule height0.3em}%
|
||||
}%
|
||||
\par\penalty10000%
|
||||
}
|
||||
|
||||
%% callable by ref
|
||||
\def\delimtitleref#1{
|
||||
\par%
|
||||
%
|
||||
\ifdefined\ianclass%
|
||||
% hyperref anchor%
|
||||
\hrefanchor%
|
||||
% define tag (for \label)%
|
||||
\xdef\tag{#1}%
|
||||
\fi%
|
||||
%
|
||||
\leavevmode%
|
||||
\raise.3em\hbox to\hsize{%
|
||||
\lower0.3em\hbox{\vrule height0.3em}%
|
||||
\hrulefill%
|
||||
\ \lower.3em\hbox{\bf #1}\ %
|
||||
\hrulefill%
|
||||
\lower0.3em\hbox{\vrule height0.3em}%
|
||||
}%
|
||||
\par\penalty10000%
|
||||
}
|
||||
|
||||
%% no title
|
||||
\def\delim{
|
||||
\par%
|
||||
\leavevmode\raise.3em\hbox to\hsize{%
|
||||
\lower0.3em\hbox{\vrule height0.3em}%
|
||||
\hrulefill%
|
||||
\lower0.3em\hbox{\vrule height0.3em}%
|
||||
}%
|
||||
\par\penalty10000%
|
||||
}
|
||||
|
||||
%% end delim
|
||||
\def\enddelim{
|
||||
\par\penalty10000%
|
||||
\leavevmode%
|
||||
\raise.3em\hbox to\hsize{%
|
||||
\vrule height0.3em\hrulefill\vrule height0.3em%
|
||||
}%
|
||||
\par%
|
||||
}
|
||||
|
||||
|
||||
%% theorem
|
||||
% counter
|
||||
\newcounter{theocount}
|
||||
% booleans (write section or subsection in equation number)
|
||||
\def\theo#1{
|
||||
\stepcounter{theocount}
|
||||
\ifdefined\ianclass
|
||||
% hyperref anchor
|
||||
\hrefanchor
|
||||
\fi
|
||||
% the number
|
||||
\def\formattheo{\thetheocount}
|
||||
% add subsection number
|
||||
\ifsubsectionintheo
|
||||
\let\tmp\formattheo
|
||||
\edef\formattheo{\thesubsectioncount.\tmp}
|
||||
\fi
|
||||
% add section number
|
||||
\ifsectionintheo
|
||||
\let\tmp\formattheo
|
||||
\edef\formattheo{\sectionprefix\thesectioncount.\tmp}
|
||||
\fi
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\formattheo}
|
||||
% write
|
||||
\delimtitle{\bf #1 \formattheo}
|
||||
}
|
||||
\let\endtheo\enddelim
|
||||
%% theorem headers with name
|
||||
\def\theoname#1#2{
|
||||
\theo{#1}\hfil({\it #2})\par\penalty10000\medskip%
|
||||
}
|
||||
|
||||
|
||||
%% qed symbol
|
||||
\def\qedsymbol{$\square$}
|
||||
\def\qed{\penalty10000\hfill\penalty10000\qedsymbol}
|
||||
|
||||
|
||||
%% compatibility with article class
|
||||
\ifdefined\ianclasstrue
|
||||
\relax
|
||||
\else
|
||||
\def\thesectioncount{\thesection}
|
||||
\def\thesubsectioncount{\thesubsection}
|
||||
\def\sectionprefix{}
|
||||
\fi
|
||||
|
||||
|
||||
%% prevent page breaks after displayed equations
|
||||
\newcount\prevpostdisplaypenalty
|
||||
\def\nopagebreakaftereq{
|
||||
\prevpostdisplaypenalty=\postdisplaypenalty
|
||||
\postdisplaypenalty=10000
|
||||
}
|
||||
%% back to previous value
|
||||
\def\restorepagebreakaftereq{
|
||||
\postdisplaypenalty=\prevpostdisplaypenalty
|
||||
}
|
||||
|
||||
|
||||
%% end
|
||||
\iantheo@defaultoptions
|
||||
\endinput
|
||||
33
doc/libs/largearray.sty
Normal file
33
doc/libs/largearray.sty
Normal file
@ -0,0 +1,33 @@
|
||||
%% Copyright 2021 Ian Jauslin
|
||||
%%
|
||||
%% Licensed under the Apache License, Version 2.0 (the "License");
|
||||
%% you may not use this file except in compliance with the License.
|
||||
%% You may obtain a copy of the License at
|
||||
%%
|
||||
%% http://www.apache.org/licenses/LICENSE-2.0
|
||||
%%
|
||||
%% Unless required by applicable law or agreed to in writing, software
|
||||
%% distributed under the License is distributed on an "AS IS" BASIS,
|
||||
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
%% See the License for the specific language governing permissions and
|
||||
%% limitations under the License.
|
||||
|
||||
%%
|
||||
%% largearray package:
|
||||
%% Array spanning the entire line
|
||||
%%
|
||||
|
||||
%% TeX format
|
||||
\NeedsTeXFormat{LaTeX2e}[1995/12/01]
|
||||
|
||||
%% package name
|
||||
\ProvidesPackage{largearray}[2016/11/10]
|
||||
|
||||
\RequirePackage{array}
|
||||
|
||||
%% array spanning the entire line
|
||||
\newlength\largearray@width
|
||||
\setlength\largearray@width\textwidth
|
||||
\addtolength\largearray@width{-10pt}
|
||||
\def\largearray{\begin{array}{@{}>{\displaystyle}l@{}}\hphantom{\hspace{\largearray@width}}\\[-.5cm]}
|
||||
\def\endlargearray{\end{array}}
|
||||
128
doc/libs/point.sty
Normal file
128
doc/libs/point.sty
Normal file
@ -0,0 +1,128 @@
|
||||
%% Copyright 2021 Ian Jauslin
|
||||
%%
|
||||
%% Licensed under the Apache License, Version 2.0 (the "License");
|
||||
%% you may not use this file except in compliance with the License.
|
||||
%% You may obtain a copy of the License at
|
||||
%%
|
||||
%% http://www.apache.org/licenses/LICENSE-2.0
|
||||
%%
|
||||
%% Unless required by applicable law or agreed to in writing, software
|
||||
%% distributed under the License is distributed on an "AS IS" BASIS,
|
||||
%% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
%% See the License for the specific language governing permissions and
|
||||
%% limitations under the License.
|
||||
|
||||
%%
|
||||
%% Points package:
|
||||
%% \point commands
|
||||
%%
|
||||
|
||||
%% TeX format
|
||||
\NeedsTeXFormat{LaTeX2e}[1995/12/01]
|
||||
|
||||
%% package name
|
||||
\ProvidesPackage{point}[2017/06/13]
|
||||
|
||||
%% options
|
||||
\newif\ifresetatsection
|
||||
\DeclareOption{reset_at_section}{\resetatsectiontrue}
|
||||
\DeclareOption{no_reset_at_section}{\resetatsectionfalse}
|
||||
\newif\ifresetatsubsection
|
||||
\DeclareOption{reset_at_subsection}{\resetatsubsectiontrue}
|
||||
\DeclareOption{no_reset_at_subsection}{\resetatsubsectionfalse}
|
||||
\newif\ifresetatsubsubsection
|
||||
\DeclareOption{reset_at_subsubsection}{\resetatsubsubsectiontrue}
|
||||
\DeclareOption{no_reset_at_subsubsection}{\resetatsubsubsectionfalse}
|
||||
\newif\ifresetattheo
|
||||
\DeclareOption{reset_at_theo}{\resetattheotrue}
|
||||
\DeclareOption{no_reset_at_theo}{\resetattheofalse}
|
||||
|
||||
\def\point@defaultoptions{
|
||||
\ExecuteOptions{reset_at_section, reset_at_subsection, reset_at_subsubsection, no_reset_at_theo}
|
||||
\ProcessOptions
|
||||
|
||||
%% reset at every new section
|
||||
\ifresetatsection
|
||||
\let\point@oldsection\section
|
||||
\gdef\section{\resetpointcounter\point@oldsection}
|
||||
\fi
|
||||
%% reset at every new subsection
|
||||
\ifresetatsubsection
|
||||
\let\point@oldsubsection\subsection
|
||||
\gdef\subsection{\resetpointcounter\point@oldsubsection}
|
||||
\fi
|
||||
%% reset at every new subsubsection
|
||||
\ifresetatsubsubsection
|
||||
\let\point@oldsubsubsection\subsubsection
|
||||
\gdef\subsubsection{\resetpointcounter\point@oldsubsubsection}
|
||||
\fi
|
||||
|
||||
%% reset at every new theorem
|
||||
\ifresetattheo
|
||||
\ifdefined\iantheotrue
|
||||
\let\point@oldtheo\theo
|
||||
\gdef\theo{\resetpointcounter\point@oldtheo}
|
||||
\fi
|
||||
\fi
|
||||
}
|
||||
|
||||
|
||||
%% point
|
||||
% counter
|
||||
\newcounter{pointcount}
|
||||
\def\point{
|
||||
\stepcounter{pointcount}
|
||||
\setcounter{subpointcount}{0}
|
||||
% hyperref anchor (only if the class is 'ian')
|
||||
\ifdefined\ifianclass
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\thepointcount}
|
||||
\fi
|
||||
% header
|
||||
\indent{\bf \thepointcount\ - }
|
||||
}
|
||||
|
||||
%% subpoint
|
||||
% counter
|
||||
\newcounter{subpointcount}
|
||||
\def\subpoint{
|
||||
\stepcounter{subpointcount}
|
||||
\setcounter{subsubpointcount}0
|
||||
% hyperref anchor (only if the class is 'ian')
|
||||
\ifdefined\ifianclass
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\thepointcount-\thesubpointcount}
|
||||
\fi
|
||||
% header
|
||||
\indent\hskip.5cm{\bf \thepointcount-\thesubpointcount\ - }
|
||||
}
|
||||
|
||||
%% subsubpoint
|
||||
% counter
|
||||
\newcounter{subsubpointcount}
|
||||
\def\subsubpoint{
|
||||
\stepcounter{subsubpointcount}
|
||||
% hyperref anchor (only if the class is 'ian')
|
||||
\ifdefined\ifianclass
|
||||
\hrefanchor
|
||||
% define tag (for \label)
|
||||
\xdef\tag{\thepointcount-\thesubpointcount-\thesubsubpointcount}
|
||||
\fi
|
||||
\indent\hskip1cm{\bf \thepointcount-\thesubpointcount-\thesubsubpointcount\ - }
|
||||
}
|
||||
|
||||
|
||||
%% reset point counters
|
||||
\def\resetpointcounter{
|
||||
\setcounter{pointcount}{0}
|
||||
\setcounter{subpointcount}{0}
|
||||
\setcounter{subsubpointcount}{0}
|
||||
}
|
||||
|
||||
|
||||
|
||||
%% end
|
||||
\point@defaultoptions
|
||||
\endinput
|
||||
3679
doc/simplesolv-doc.tex
Normal file
3679
doc/simplesolv-doc.tex
Normal file
File diff suppressed because it is too large
Load Diff
949
src/anyeq.jl
Normal file
949
src/anyeq.jl
Normal file
@ -0,0 +1,949 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# interpolation
|
||||
@everywhere mutable struct Anyeq_approx
|
||||
aK::Float64
|
||||
bK::Float64
|
||||
gK::Float64
|
||||
aL1::Float64
|
||||
bL1::Float64
|
||||
aL2::Float64
|
||||
bL2::Float64
|
||||
gL2::Float64
|
||||
aL3::Float64
|
||||
bL3::Float64
|
||||
gL3::Float64
|
||||
end
|
||||
|
||||
# compute energy for a given rho
|
||||
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
|
||||
function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
|
||||
@printf("% .15e % .15e\n",E,error)
|
||||
end
|
||||
|
||||
# compute energy as a function of rho
|
||||
function anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
|
||||
# save result from each task
|
||||
es=Array{Float64,1}(undef,length(rhos))
|
||||
err=Array{Float64,1}(undef,length(rhos))
|
||||
|
||||
## spawn workers
|
||||
# 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]=zeros(0)
|
||||
end
|
||||
for j in 1:length(rhos)
|
||||
append!(work[(j-1)%nw+1],j)
|
||||
end
|
||||
|
||||
count=0
|
||||
# for each worker
|
||||
@sync for p in 1:nw
|
||||
# for each task
|
||||
@async for j in work[p]
|
||||
count=count+1
|
||||
if count>=nw
|
||||
progress(count,length(rhos),10000)
|
||||
end
|
||||
# run the task
|
||||
(u,es[j],err[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
end
|
||||
end
|
||||
|
||||
for j in 1:length(rhos)
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j])
|
||||
end
|
||||
end
|
||||
|
||||
# compute energy as a function of rho
|
||||
# initialize with previous rho
|
||||
function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
u0s=anyeq_init_medeq([rhos[1]],N,J,k,a0,v,maxiter,tolerance)
|
||||
u=u0s[1]
|
||||
|
||||
for j in 1:length(rhos)
|
||||
progress(j,length(rhos),10000)
|
||||
# run the task
|
||||
# init Newton from previous rho
|
||||
(u,E,error)=anyeq_hatu(u,P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],E,error)
|
||||
# abort when the error gets too big
|
||||
if error>tolerance
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
# compute energy as a function of rho
|
||||
# initialize with next rho
|
||||
function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
u0s=anyeq_init_medeq([rhos[length(rhos)]],N,J,k,a0,v,maxiter,tolerance)
|
||||
u=u0s[1]
|
||||
|
||||
for j in 1:length(rhos)
|
||||
progress(j,length(rhos),10000)
|
||||
# run the task
|
||||
# init Newton from previous rho
|
||||
(u,E,error)=anyeq_hatu(u,P,N,J,rhos[length(rhos)+1-j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
@printf("% .15e % .15e % .15e\n",rhos[length(rhos)+1-j],real(E),error)
|
||||
# abort when the error gets too big
|
||||
if error>tolerance
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# compute u(k)
|
||||
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
|
||||
function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# init vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
# order k's in increasing order
|
||||
@printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],u[(J-1-zeta)*N+j])
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# compute u(x)
|
||||
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
|
||||
function anyeq_ux(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile)
|
||||
# init vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
ux=0.
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
ux+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*u[zeta*N+j]*sin(k[zeta*N+j]*x)
|
||||
end
|
||||
end
|
||||
@printf("% .15e % .15e % .15e\n",x,real(ux),imag(ux))
|
||||
end
|
||||
end
|
||||
|
||||
# compute condensate fraction for a given rho
|
||||
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
|
||||
function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
|
||||
# compute eta
|
||||
eta=anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
|
||||
|
||||
@printf("% .15e % .15e\n",eta,error)
|
||||
end
|
||||
|
||||
# condensate fraction as a function of rho
|
||||
function anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
## spawn workers
|
||||
# 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]=zeros(0)
|
||||
end
|
||||
for j in 1:length(rhos)
|
||||
append!(work[(j-1)%nw+1],j)
|
||||
end
|
||||
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
|
||||
# compute u
|
||||
us=Array{Array{Float64,1}}(undef,length(rhos))
|
||||
errs=Array{Float64,1}(undef,length(rhos))
|
||||
count=0
|
||||
# for each worker
|
||||
@sync for p in 1:nw
|
||||
# for each task
|
||||
@async for j in work[p]
|
||||
count=count+1
|
||||
if count>=nw
|
||||
progress(count,length(rhos),10000)
|
||||
end
|
||||
# run the task
|
||||
(us[j],E,errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
end
|
||||
end
|
||||
|
||||
# compute eta
|
||||
etas=Array{Float64}(undef,length(rhos))
|
||||
count=0
|
||||
# for each worker
|
||||
@sync for p in 1:nw
|
||||
# for each task
|
||||
@async for j in work[p]
|
||||
count=count+1
|
||||
if count>=nw
|
||||
progress(count,length(rhos),10000)
|
||||
end
|
||||
# run the task
|
||||
etas[j]=anyeq_eta(us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
|
||||
end
|
||||
end
|
||||
|
||||
for j in 1:length(rhos)
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j])
|
||||
end
|
||||
end
|
||||
|
||||
# compute the momentum distribution for a given rho
|
||||
# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm)
|
||||
function anyeq_momentum_distribution(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
|
||||
# compute M
|
||||
M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
|
||||
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
# order k's in increasing order
|
||||
@printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],M[(J-1-zeta)*N+j])
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# 2 point correlation function
|
||||
function anyeq_2pt_correlation(minlrho,nlrho,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile)
|
||||
# init vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
# init Abar
|
||||
if savefile!=""
|
||||
Abar=anyeq_read_Abar(savefile,P,N,J)
|
||||
else
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
# compute u and some useful integrals
|
||||
(u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
(S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
|
||||
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G)
|
||||
@printf("% .15e % .15e\n",x,C2)
|
||||
end
|
||||
end
|
||||
|
||||
# compute Abar
|
||||
function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
if approx.bL3==0.
|
||||
return []
|
||||
end
|
||||
|
||||
out=Array{Array{Float64,5}}(undef,J*N)
|
||||
|
||||
## spawn workers
|
||||
# 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]=zeros(0)
|
||||
end
|
||||
for j in 1:N*J
|
||||
append!(work[(j-1)%nw+1],j)
|
||||
end
|
||||
|
||||
count=0
|
||||
# for each worker
|
||||
@sync for p in 1:nw
|
||||
# for each task
|
||||
@async for j in work[p]
|
||||
count=count+1
|
||||
if count>=nw
|
||||
progress(count,N*J,10000)
|
||||
end
|
||||
# run the task
|
||||
out[j]=remotecall_fetch(barAmat,workers()[p],k[j],weights,taus,T,P,N,J,2,2,2,2,2)
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# initialize computation
|
||||
@everywhere function anyeq_init(taus,P,N,J,v)
|
||||
# Gauss-Legendre weights
|
||||
weights=gausslegendre(N)
|
||||
|
||||
# initialize vectors V,k
|
||||
V=Array{Float64}(undef,J*N)
|
||||
k=Array{Float64}(undef,J*N)
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
xj=weights[1][j]
|
||||
# set kj
|
||||
k[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1])
|
||||
# set v
|
||||
V[zeta*N+j]=v(k[zeta*N+j])
|
||||
end
|
||||
end
|
||||
# potential at 0
|
||||
V0=v(0)
|
||||
|
||||
# initialize matrix A
|
||||
T=chebyshev_polynomials(P)
|
||||
A=Amat(k,weights,taus,T,P,N,J,2,2)
|
||||
|
||||
# compute Upsilon
|
||||
# Upsilonmat does not use splines, so increase precision
|
||||
weights_plus=gausslegendre(N*J)
|
||||
Upsilon=Upsilonmat(k,v,weights_plus)
|
||||
Upsilon0=Upsilon0mat(k,v,weights_plus)
|
||||
|
||||
return(weights,T,k,V,V0,A,Upsilon,Upsilon0)
|
||||
end
|
||||
|
||||
# compute initial guess from medeq
|
||||
@everywhere function anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
us_medeq=Array{Array{Float64,1}}(undef,length(rhos))
|
||||
u0s=Array{Array{Float64,1}}(undef,length(rhos))
|
||||
|
||||
weights_medeq=gausslegendre(N*J)
|
||||
|
||||
(us_medeq[1],E,err)=easyeq_hatu(easyeq_init_u(a0,J*N,weights_medeq),J*N,rhos[1],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.))
|
||||
u0s[1]=easyeq_to_anyeq(us_medeq[1],weights_medeq,k,N,J)
|
||||
if err>tolerance
|
||||
print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[1],"\n")
|
||||
end
|
||||
|
||||
for j in 2:length(rhos)
|
||||
(us_medeq[j],E,err)=easyeq_hatu(us_medeq[j-1],J*N,rhos[j],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.))
|
||||
u0s[j]=easyeq_to_anyeq(us_medeq[j],weights_medeq,k,N,J)
|
||||
|
||||
if err>tolerance
|
||||
print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n")
|
||||
end
|
||||
end
|
||||
|
||||
return u0s
|
||||
end
|
||||
|
||||
# interpolate the solution of medeq to an input for anyeq
|
||||
@everywhere function easyeq_to_anyeq(u_simple,weights,k,N,J)
|
||||
# reorder u_simple, which is evaluated at (1-x_j)/(1+x_j) with x_j\in[-1,1]
|
||||
u_s=zeros(Float64,length(u_simple))
|
||||
k_s=Array{Float64}(undef,length(u_simple))
|
||||
for j in 1:length(u_simple)
|
||||
xj=weights[1][j]
|
||||
k_s[length(u_simple)-j+1]=(1-xj)/(1+xj)
|
||||
u_s[length(u_simple)-j+1]=u_simple[j]
|
||||
end
|
||||
|
||||
# initialize U
|
||||
u=zeros(Float64,J*N)
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
u[zeta*N+j]=linear_interpolation(k[zeta*N+j],k_s,u_s)
|
||||
end
|
||||
end
|
||||
|
||||
return u
|
||||
end
|
||||
|
||||
|
||||
# compute u using chebyshev expansions
|
||||
@everywhere function anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx)
|
||||
# init
|
||||
# rescale by rho (that's how u is defined)
|
||||
u=rho*u0
|
||||
|
||||
# quantify relative error
|
||||
error=-1.
|
||||
|
||||
# run Newton algorithm
|
||||
for i in 1:maxiter-1
|
||||
(S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
|
||||
new=u-inv(anyeq_DXi(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y)
|
||||
|
||||
error=norm(new-u)/norm(u)
|
||||
if(error<tolerance)
|
||||
u=new
|
||||
break
|
||||
end
|
||||
|
||||
u=new
|
||||
end
|
||||
|
||||
energy=rho/2*V0-avg_v_chebyshev(u,Upsilon0,k,taus,weights,N,J)/2
|
||||
return(u/rho,energy,error)
|
||||
end
|
||||
|
||||
|
||||
# save Abar
|
||||
function anyeq_save_Abar(taus,P,N,J,v,approx)
|
||||
# initialize vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# init Abar
|
||||
Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx)
|
||||
|
||||
# print params
|
||||
@printf("## P=%d N=%d J=%d\n",P,N,J)
|
||||
|
||||
for i in 1:N*J
|
||||
for j1 in 1:(P+1)*J
|
||||
for j2 in 1:(P+1)*J
|
||||
for j3 in 1:(P+1)*J
|
||||
for j4 in 1:(P+1)*J
|
||||
for j5 in 1:(P+1)*J
|
||||
@printf("% .15e\n",Abar[i][j1,j2,j3,j4,j5])
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
# read Abar
|
||||
function anyeq_read_Abar(savefile,P,N,J)
|
||||
# open file
|
||||
ff=open(savefile,"r")
|
||||
# read all lines
|
||||
lines=readlines(ff)
|
||||
close(ff)
|
||||
|
||||
# init
|
||||
Abar=Array{Array{Float64,5}}(undef,N*J)
|
||||
for i in 1:N*J
|
||||
Abar[i]=Array{Float64,5}(undef,(P+1)*J,(P+1)*J,(P+1)*J,(P+1)*J,(P+1)*J)
|
||||
end
|
||||
|
||||
i=1
|
||||
j1=1
|
||||
j2=1
|
||||
j3=1
|
||||
j4=1
|
||||
j5=0
|
||||
for l in 1:length(lines)
|
||||
# drop comments
|
||||
if lines[l]!="" && lines[l][1]!='#'
|
||||
# increment counters
|
||||
if j5<(P+1)*J
|
||||
j5+=1
|
||||
else
|
||||
j5=1
|
||||
if j4<(P+1)*J
|
||||
j4+=1
|
||||
else
|
||||
j4=1
|
||||
if j3<(P+1)*J
|
||||
j3+=1
|
||||
else
|
||||
j3=1
|
||||
if j2<(P+1)*J
|
||||
j2+=1
|
||||
else
|
||||
j2=1
|
||||
if j1<(P+1)*J
|
||||
j1+=1
|
||||
else
|
||||
j1=1
|
||||
if i<N*J
|
||||
i+=1
|
||||
else
|
||||
print(stderr,"error: too many lines in savefile\n")
|
||||
exit()
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
Abar[i][j1,j2,j3,j4,j5]=parse(Float64,lines[l])
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
return Abar
|
||||
|
||||
end
|
||||
|
||||
|
||||
|
||||
# Xi
|
||||
# takes the vector of kj's and xn's as input
|
||||
@everywhere function anyeq_Xi(U,X,Y)
|
||||
return U-(Y.+1)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))
|
||||
end
|
||||
|
||||
# DXi
|
||||
# takes the vector of kj's as input
|
||||
@everywhere function anyeq_DXi(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
|
||||
out=Array{Float64,2}(undef,N*J,N*J)
|
||||
for zetapp in 0:J-1
|
||||
for n in 1:N
|
||||
one=zeros(Int64,J*N)
|
||||
one[zetapp*N+n]=1
|
||||
|
||||
# Chebyshev expansion of U
|
||||
FU=chebyshev(U,taus,weights,P,N,J,2)
|
||||
|
||||
dS=-conv_one_v_chebyshev(n,zetapp,Upsilon,k,taus,weights,N,J)/rho
|
||||
|
||||
dE=-avg_one_v_chebyshev(n,zetapp,Upsilon0,k,taus,weights,N)/rho
|
||||
|
||||
UU=conv_chebyshev(FU,FU,A)
|
||||
|
||||
dII=zeros(Float64,N*J)
|
||||
if approx.gL2!=0.
|
||||
if approx.bL2!=0.
|
||||
dII+=approx.gL2*approx.bL2*(conv_one_chebyshev(n,zetapp,chebyshev(U.*S,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(one.*S+U.*dS,taus,weights,P,N,J,2),A)/rho)
|
||||
end
|
||||
if approx.bL2!=1.
|
||||
dII+=approx.gL2*(1-approx.bL2)*(dE/rho*UU+2*E/rho*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2))
|
||||
end
|
||||
end
|
||||
|
||||
dJJ=zeros(Float64,J*N)
|
||||
if approx.gL3!=0.
|
||||
if approx.bL3!=0.
|
||||
FS=chebyshev(S,taus,weights,P,N,J,2)
|
||||
dFU=chebyshev(one,taus,weights,P,N,J,2)
|
||||
dFS=chebyshev(dS,taus,weights,P,N,J,2)
|
||||
dJJ+=approx.gL3*approx.bL3*(4*double_conv_S_chebyshev(FU,FU,FU,dFU,FS,Abar)+double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar))
|
||||
end
|
||||
if approx.bL3!=1.
|
||||
dJJ+=approx.gL3*(1-approx.bL3)*(dE*(UU/rho).^2+4*E*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2).*UU/rho^2)
|
||||
end
|
||||
end
|
||||
|
||||
dG=zeros(Float64,N*J)
|
||||
if approx.aK!=0. && approx.gK!=0.
|
||||
if approx.bK!=0.
|
||||
dG+=approx.aK*approx.gK*approx.bK*(conv_one_chebyshev(n,zetapp,chebyshev(2*S.*U,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*S.*one+2*dS.*U,taus,weights,P,N,J,2),A)/rho)
|
||||
end
|
||||
if approx.bK!=1.
|
||||
dG+=approx.aK*approx.gK*(1-approx.bK)*(2*dE*UU/rho+4*E*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2)/rho)
|
||||
end
|
||||
end
|
||||
if approx.aL1!=0.
|
||||
if approx.bL1!=0.
|
||||
dG-=approx.aL1*approx.bL1*(conv_one_chebyshev(n,zetapp,chebyshev(S.*(U.^2),taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*S.*U.*one+dS.*(U.^2),taus,weights,P,N,J,2),A)/rho)
|
||||
end
|
||||
if approx.bL1!=1.
|
||||
dG-=approx.aL1*(1-approx.bL1)*(E/rho*conv_one_chebyshev(n,zetapp,chebyshev((U.^2),taus,weights,P,N,J,2),taus,A,weights,P,N,J,2)+conv_chebyshev(FU,chebyshev(2*E*U.*one+dE*(U.^2),taus,weights,P,N,J,2),A)/rho)
|
||||
end
|
||||
end
|
||||
if approx.aL2!=0. && approx.gL2!=0.
|
||||
dG+=approx.aL2*(conv_one_chebyshev(n,zetapp,chebyshev(2*II.*U,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*dII.*U+2*II.*one,taus,weights,P,N,J,2),A)/rho)
|
||||
end
|
||||
if approx.aL3!=0. && approx.gL3!=0.
|
||||
dG-=approx.aL3*(conv_one_chebyshev(n,zetapp,chebyshev(JJ/2,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho)
|
||||
end
|
||||
|
||||
dsK=zeros(Float64,N*J)
|
||||
if approx.gK!=0.
|
||||
if approx.bK!=0.
|
||||
dsK+=approx.gK*approx.bK*dS
|
||||
end
|
||||
if approx.bK!=1.
|
||||
dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J)
|
||||
end
|
||||
else
|
||||
end
|
||||
dsL1=zeros(Float64,N*J)
|
||||
if approx.bL1!=0.
|
||||
dsL1+=approx.bL1*dS
|
||||
end
|
||||
if approx.bL1!=1.
|
||||
dsL1+=(1-approx.bL1)*dE*ones(Float64,N*J)
|
||||
end
|
||||
|
||||
dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1
|
||||
dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1
|
||||
|
||||
out[:,zetapp*N+n]=one+((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# compute S,E,I,J,X and Y
|
||||
@everywhere function anyeq_SEIJGXY(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
|
||||
# Chebyshev expansion of U
|
||||
FU=chebyshev(U,taus,weights,P,N,J,2)
|
||||
|
||||
S=v-conv_v_chebyshev(U,Upsilon,k,taus,weights,N,J)/rho
|
||||
E=v0-avg_v_chebyshev(U,Upsilon0,k,taus,weights,N,J)/rho
|
||||
|
||||
UU=conv_chebyshev(FU,FU,A)
|
||||
|
||||
II=zeros(Float64,N*J)
|
||||
if approx.gL2!=0.
|
||||
if approx.bL2!=0.
|
||||
II+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*S,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.bL2!=1.
|
||||
II+=approx.gL2*(1-approx,bL2)*E/rho*UU
|
||||
end
|
||||
end
|
||||
|
||||
JJ=zeros(Float64,N*J)
|
||||
if approx.gL3!=0.
|
||||
if approx.bL3!=0.
|
||||
FS=chebyshev(S,taus,weights,P,N,J,2)
|
||||
JJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,FS,Abar)
|
||||
end
|
||||
if approx.bL3!=1.
|
||||
JJ+=approx.gL3*(1-approx.bL3)*E*(UU/rho).^2
|
||||
end
|
||||
end
|
||||
|
||||
G=zeros(Float64,N*J)
|
||||
if approx.aK!=0. && approx.gK!=0.
|
||||
if approx.bK!=0.
|
||||
G+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*S.*U,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.bK!=1.
|
||||
G+=approx.aK*approx.gK*(1-approx.bK)*2*E*UU/rho
|
||||
end
|
||||
end
|
||||
if approx.aL1!=0.
|
||||
if approx.bL1!=0.
|
||||
G-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(S.*(U.^2),taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.bL1!=1.
|
||||
G-=approx.aL1*(1-approx.bL1)*E/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A)
|
||||
end
|
||||
end
|
||||
if approx.aL2!=0. && approx.gL2!=0.
|
||||
G+=approx.aL2*conv_chebyshev(FU,chebyshev(2*II.*U,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.aL3!=0 && approx.gL3!=0.
|
||||
G-=approx.aL3*conv_chebyshev(FU,chebyshev(JJ/2,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
|
||||
sK=zeros(Float64,N*J)
|
||||
if approx.gK!=0.
|
||||
if approx.bK!=0.
|
||||
sK+=approx.gK*approx.bK*S
|
||||
end
|
||||
if approx.bK!=1.
|
||||
sK+=approx.gK*(1-approx.bK)*E*ones(Float64,N*J)
|
||||
end
|
||||
end
|
||||
|
||||
sL1=zeros(Float64,N*J)
|
||||
if approx.bL1!=0.
|
||||
sL1+=approx.bL1*S
|
||||
end
|
||||
if approx.bL1!=1.
|
||||
sL1+=(1-approx.bL1)*E*ones(Float64,N*J)
|
||||
end
|
||||
|
||||
X=(k.^2/2+rho*(sK-sL1+II))./sL1/rho
|
||||
Y=(S-sL1+G+JJ/2)./sL1
|
||||
|
||||
return(S,E,II,JJ,X,Y,sL1,sK,G)
|
||||
end
|
||||
|
||||
# condensate fraction
|
||||
@everywhere function anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
|
||||
# compute dXi/dmu
|
||||
(S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
|
||||
dXidmu=(Y.+1)./(rho*sL1)./(2*(X.+1).^2).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1).^2 ./((X.+1).^4)./(rho*sL1).*dotdPhi((Y.+1)./(X.+1).^2)
|
||||
|
||||
# compute eta
|
||||
du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu
|
||||
eta=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2
|
||||
|
||||
return eta
|
||||
end
|
||||
|
||||
# momentum distribution
|
||||
@everywhere function anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx)
|
||||
# compute dXi/dlambda (without delta functions)
|
||||
(S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx)
|
||||
dXidlambda=-(2*pi)^3*2*u./sL1.*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2))
|
||||
|
||||
# approximation for delta function (without Kronecker deltas)
|
||||
delta=Array{Float64}(undef,J*N)
|
||||
for zeta in 0:J-1
|
||||
for n in 1:N
|
||||
delta[zeta*N+n]=2/pi^2/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2*k[zeta*N+n]^2)
|
||||
end
|
||||
end
|
||||
|
||||
# compute dXidu
|
||||
dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))
|
||||
|
||||
M=Array{Float64}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
# du/dlambda
|
||||
du=dXidu[:,i]*dXidlambda[i]*delta[i]
|
||||
|
||||
# compute M
|
||||
M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2
|
||||
end
|
||||
|
||||
return M
|
||||
end
|
||||
|
||||
|
||||
# correlation function
|
||||
@everywhere function anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G)
|
||||
# initialize dV
|
||||
dV=Array{Float64}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
if x>0
|
||||
dV[i]=sin(k[i]*x)/(k[i]*x)*hann(k[i],windowL)
|
||||
else
|
||||
dV[i]=hann(k[i],windowL)
|
||||
end
|
||||
end
|
||||
dV0=1.
|
||||
|
||||
# compute dUpsilon
|
||||
# Upsilonmat does not use splines, so increase precision
|
||||
weights_plus=gausslegendre(N*J)
|
||||
dUpsilon=Upsilonmat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus)
|
||||
dUpsilon0=Upsilon0mat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus)
|
||||
|
||||
du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(x,rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
|
||||
# rescale rho
|
||||
du=du/rho
|
||||
|
||||
C2=rho^2*(1-integrate_f_chebyshev(s->1.,u.*dV+V.*du,k,taus,weights,N,J))
|
||||
|
||||
return C2
|
||||
end
|
||||
|
||||
# derivative of Xi with respect to v in the direction sin(kx)/kx
|
||||
@everywhere function anyeq_dXidv(x,U,rho,k,taus,dv,dv0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)
|
||||
# Chebyshev expansion of U
|
||||
FU=chebyshev(U,taus,weights,P,N,J,2)
|
||||
|
||||
dS=dv-conv_v_chebyshev(U,dUpsilon,k,taus,weights,N,J)/rho
|
||||
dE=dv0-avg_v_chebyshev(U,dUpsilon0,k,taus,weights,N,J)/rho
|
||||
|
||||
UU=conv_chebyshev(FU,FU,A)
|
||||
|
||||
dII=zeros(Float64,N*J)
|
||||
if approx.gL2!=0.
|
||||
if approx.bL2!=0.
|
||||
dII+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*dS,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.bL2!=1.
|
||||
dII+=approx.gL2*(1-approx,bL2)*dE/rho*UU
|
||||
end
|
||||
end
|
||||
|
||||
dJJ=zeros(Float64,J*N)
|
||||
if approx.gL3!=0.
|
||||
if approx.bL3!=0.
|
||||
dFS=chebyshev(dS,taus,weights,P,N,J,2)
|
||||
dJJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar)
|
||||
end
|
||||
if approx.bL3!=1.
|
||||
dJJ=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2
|
||||
end
|
||||
end
|
||||
|
||||
dG=zeros(Float64,N*J)
|
||||
if approx.aK!=0. && approx.gK!=0.
|
||||
if approx.bK!=0.
|
||||
dG+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*dS.*U,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.bK!=1.
|
||||
dG+=approx.aK*approx.gK*(1-approx.bK)*2*dE*UU/rho
|
||||
end
|
||||
end
|
||||
if approx.aL1!=0.
|
||||
if approx.bL1!=0.
|
||||
dG-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(dS.*(U.^2),taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.bL1!=1.
|
||||
dG-=approx.aL1*(1-approx.bL1)*dE/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A)
|
||||
end
|
||||
end
|
||||
if approx.aL2!=0. && approx.gL2!=0.
|
||||
dG+=approx.aL2*conv_chebyshev(FU,chebyshev(2*dII.*U,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
if approx.aL3!=0. && approx.gL3!=0.
|
||||
dG-=approx.aL3*conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho
|
||||
end
|
||||
|
||||
dsK=zeros(Float64,N*J)
|
||||
if approx.gK!=0.
|
||||
if approx.bK!=0.
|
||||
dsK+=approx.gK*approx.bK*dS
|
||||
end
|
||||
if approx.bK!=1.
|
||||
dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J)
|
||||
end
|
||||
end
|
||||
dsL1=zeros(Float64,N*J)
|
||||
if approx.bL1!=0.
|
||||
dsL1+=approx.bL1*dS
|
||||
end
|
||||
if approx.bL1!=1.
|
||||
dsL1+=(1-approx.bL1)*dE*ones(N*J)
|
||||
end
|
||||
|
||||
dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1
|
||||
dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1
|
||||
|
||||
out=((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2)
|
||||
return out
|
||||
end
|
||||
279
src/chebyshev.jl
Normal file
279
src/chebyshev.jl
Normal file
@ -0,0 +1,279 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# Chebyshev expansion
|
||||
@everywhere function chebyshev(a,taus,weights,P,N,J,nu)
|
||||
out=zeros(Float64,J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for j in 1:N
|
||||
out[zeta*(P+1)+n+1]+=(2-(n==0 ? 1 : 0))/2*weights[2][j]*cos(n*pi*(1+weights[1][j])/2)*a[zeta*N+j]/(1-(taus[zeta+2]-taus[zeta+1])/2*sin(pi*weights[1][j]/2)+(taus[zeta+2]+taus[zeta+1])/2)^nu
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# evaluate function from Chebyshev expansion
|
||||
@everywhere function chebyshev_eval(Fa,x,taus,chebyshev,P,J,nu)
|
||||
# change variable
|
||||
tau=(1-x)/(1+x)
|
||||
|
||||
out=0.
|
||||
for zeta in 0:J-1
|
||||
# check that the spline is right
|
||||
if tau<taus[zeta+2] && tau>=taus[zeta+1]
|
||||
for n in 0:P
|
||||
out+=Fa[zeta*(P+1)+n+1]*chebyshev[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))
|
||||
end
|
||||
break
|
||||
end
|
||||
end
|
||||
|
||||
return (1+tau)^nu*out
|
||||
end
|
||||
|
||||
|
||||
# convolution
|
||||
# input the Chebyshev expansion of a and b, as well as the A matrix
|
||||
@everywhere function conv_chebyshev(Fa,Fb,A)
|
||||
out=zeros(Float64,length(A))
|
||||
for i in 1:length(A)
|
||||
out[i]=dot(Fa,A[i]*Fb)
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# <ab>
|
||||
@everywhere function avg_chebyshev(Fa,Fb,A0)
|
||||
return dot(Fa,A0*Fb)
|
||||
end
|
||||
|
||||
# 1_n * a
|
||||
@everywhere function conv_one_chebyshev(n,zetapp,Fa,A,taus,weights,P,N,J,nu1)
|
||||
out=zeros(Float64,N*J)
|
||||
for m in 1:N*J
|
||||
for l in 0:P
|
||||
for p in 1:J*(P+1)
|
||||
out[m]+=(2-(l==0 ? 1 : 0))/2*weights[2][n]*cos(l*pi*(1+weights[1][n])/2)/(1-(taus[zetapp+2]-taus[zetapp+1])/2*sin(pi*weights[1][n]/2)+(taus[zetapp+2]+taus[zetapp+1])/2)^nu1*A[m][zetapp*(P+1)+l+1,p]*Fa[p]
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
# a * v
|
||||
@everywhere function conv_v_chebyshev(a,Upsilon,k,taus,weights,N,J)
|
||||
out=zeros(Float64,J*N)
|
||||
for i in 1:J*N
|
||||
for zetap in 0:J-1
|
||||
for j in 1:N
|
||||
xj=weights[1][j]
|
||||
out[i]+=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][j]*cos(pi*xj/2)*(1+k[zetap*N+j])^2*k[zetap*N+j]*a[zetap*N+j]*Upsilon[zetap*N+j][i]
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function conv_one_v_chebyshev(n,zetap,Upsilon,k,taus,weights,N,J)
|
||||
out=zeros(Float64,J*N)
|
||||
xj=weights[1][n]
|
||||
for i in 1:J*N
|
||||
out[i]=(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]*Upsilon[zetap*N+n][i]
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# <av>
|
||||
@everywhere function avg_v_chebyshev(a,Upsilon0,k,taus,weights,N,J)
|
||||
out=0.
|
||||
for zetap in 0:J-1
|
||||
for j in 1:N
|
||||
xj=weights[1][j]
|
||||
out+=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][j]*cos(pi*xj/2)*(1+k[zetap*N+j])^2*k[zetap*N+j]*a[zetap*N+j]*Upsilon0[zetap*N+j]
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
# <1_nv>
|
||||
@everywhere function avg_one_v_chebyshev(n,zetap,Upsilon0,k,taus,weights,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]
|
||||
end
|
||||
|
||||
# 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)
|
||||
out=zeros(Float64,length(Abar))
|
||||
for i in 1:length(Abar)
|
||||
for j1 in 1:length(FU1)
|
||||
for j2 in 1:length(FU2)
|
||||
for j3 in 1:length(FU3)
|
||||
for j4 in 1:length(FU4)
|
||||
for j5 in 1:length(FU5)
|
||||
out[i]+=Abar[i][j1,j2,j3,j4,j5]*FU1[j1]*FU2[j2]*FU3[j3]*FU4[j4]*FU5[j5]
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
|
||||
# compute A
|
||||
@everywhere function Amat(k,weights,taus,T,P,N,J,nua,nub)
|
||||
out=Array{Array{Float64,2},1}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=zeros(Float64,J*(P+1),J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
for m in 0:P
|
||||
out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=1/(pi^2*k[i])*integrate_legendre(tau->(1-tau)/(1+tau)^(3-nua)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))*(alpham(k[i],tau)>taus[zetap+2] || alphap(k[i],tau)<taus[zetap+1] ? 0. : integrate_legendre(sigma->(1-sigma)/(1+sigma)^(3-nub)*T[m+1]((2*sigma-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1],alpham(k[i],tau)),min(taus[zetap+2],alphap(k[i],tau)),weights)),taus[zeta+1],taus[zeta+2],weights)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# compute Upsilon
|
||||
@everywhere function Upsilonmat(k,v,weights)
|
||||
out=Array{Array{Float64,1},1}(undef,length(k))
|
||||
for i in 1:length(k)
|
||||
out[i]=Array{Float64,1}(undef,length(k))
|
||||
for j in 1:length(k)
|
||||
out[i][j]=integrate_legendre(s->s*v(s)/k[j],abs(k[j]-k[i]),k[j]+k[i],weights)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function Upsilon0mat(k,v,weights)
|
||||
out=Array{Float64,1}(undef,length(k))
|
||||
for j in 1:length(k)
|
||||
out[j]=2*k[j]*v(k[j])
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# alpha_-
|
||||
@everywhere function alpham(k,t)
|
||||
return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t))
|
||||
end
|
||||
# alpha_+
|
||||
@everywhere function alphap(k,t)
|
||||
return (1-abs(k-(1-t)/(1+t)))/(1+abs(k-(1-t)/(1+t)))
|
||||
end
|
||||
|
||||
|
||||
# compute \bar A
|
||||
@everywhere function barAmat(k,weights,taus,T,P,N,J,nu1,nu2,nu3,nu4,nu5)
|
||||
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 n1 in 0:P
|
||||
for zeta2 in 0:J-1
|
||||
for n2 in 0:P
|
||||
for zeta3 in 0:J-1
|
||||
for n3 in 0:P
|
||||
for zeta4 in 0:J-1
|
||||
for n4 in 0:P
|
||||
for zeta5 in 0:J-1
|
||||
for n5 in 0:P
|
||||
out[zeta1*(P+1)+n1+1,
|
||||
zeta2*(P+1)+n2+1,
|
||||
zeta3*(P+1)+n3+1,
|
||||
zeta4*(P+1)+n4+1,
|
||||
zeta5*(P+1)+n5+1]=1/((2*pi)^5*k^2)*integrate_legendre(tau->barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5),taus[zeta1+1],taus[zeta1+2],weights)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
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)
|
||||
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)
|
||||
else
|
||||
return 0.
|
||||
end
|
||||
end
|
||||
@everywhere function barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5)
|
||||
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
|
||||
@everywhere function barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5)
|
||||
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)
|
||||
else
|
||||
return 0.
|
||||
end
|
||||
end
|
||||
@everywhere function barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5)
|
||||
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
|
||||
@everywhere function barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5)
|
||||
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])
|
||||
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]))
|
||||
else
|
||||
return 0.
|
||||
end
|
||||
end
|
||||
# R(s,t,s',t,theta,k)
|
||||
@everywhere function barAmat_R(s,t,sp,tp,theta,k)
|
||||
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
|
||||
|
||||
# compute Chebyshev polynomials
|
||||
@everywhere function chebyshev_polynomials(P)
|
||||
T=Array{Polynomial}(undef,P+1)
|
||||
T[1]=Polynomial([1])
|
||||
T[2]=Polynomial([0,1])
|
||||
for n in 1:P-1
|
||||
# T_n
|
||||
T[n+2]=2*T[2]*T[n+1]-T[n]
|
||||
end
|
||||
|
||||
return T
|
||||
end
|
||||
|
||||
# compute \int f*u dk/(2*pi)^3
|
||||
@everywhere function integrate_f_chebyshev(f,u,k,taus,weights,N,J)
|
||||
out=0.
|
||||
for zeta in 0:J-1
|
||||
for i in 1:N
|
||||
out+=(taus[zeta+2]-taus[zeta+1])/(16*pi)*weights[2][i]*cos(pi*weights[1][i]/2)*(1+k[zeta*N+i])^2*k[zeta*N+i]^2*u[zeta*N+i]*f(k[zeta*N+i])
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
@everywhere function inverse_fourier_chebyshev(u,x,k,taus,weights,N,J)
|
||||
out=0.
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
out+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*u[zeta*N+j]*sin(k[zeta*N+j]*x)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
411
src/easyeq.jl
Normal file
411
src/easyeq.jl
Normal file
@ -0,0 +1,411 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# interpolation
|
||||
@everywhere mutable struct Easyeq_approx
|
||||
bK::Float64
|
||||
bL::Float64
|
||||
end
|
||||
|
||||
# compute energy
|
||||
function easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx)
|
||||
# compute gaussian quadrature weights
|
||||
weights=gausslegendre(order)
|
||||
|
||||
# compute initial guess from previous rho
|
||||
(u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
|
||||
for j in 2:nlrho_init
|
||||
rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
|
||||
(u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
|
||||
end
|
||||
|
||||
# print energy
|
||||
@printf("% .15e % .15e\n",real(E),err)
|
||||
end
|
||||
|
||||
# compute energy as a function of rho
|
||||
function easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,approx)
|
||||
# compute gaussian quadrature weights
|
||||
weights=gausslegendre(order)
|
||||
# init u
|
||||
u=easyeq_init_u(a0,order,weights)
|
||||
|
||||
for j in 1:length(rhos)
|
||||
# compute u (init newton with previously computed u)
|
||||
(u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx)
|
||||
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],real(E),err)
|
||||
|
||||
end
|
||||
end
|
||||
|
||||
# compute u(k)
|
||||
function easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx)
|
||||
weights=gausslegendre(order)
|
||||
|
||||
# compute initial guess from previous rho
|
||||
(u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
|
||||
for j in 2:nlrho_init
|
||||
rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
|
||||
(u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
|
||||
end
|
||||
|
||||
for i in 1:order
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
@printf("% .15e % .15e\n",k,real(u[i]))
|
||||
end
|
||||
end
|
||||
|
||||
# compute u(x)
|
||||
function easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx)
|
||||
weights=gausslegendre(order)
|
||||
|
||||
# compute initial guess from previous rho
|
||||
(u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
|
||||
for j in 2:nlrho_init
|
||||
rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
|
||||
(u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
|
||||
end
|
||||
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
@printf("% .15e % .15e\n",x,real(easyeq_u_x(x,u,weights)))
|
||||
end
|
||||
end
|
||||
|
||||
# compute 2u(x)-rho u*u(x)
|
||||
function easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx)
|
||||
weights=gausslegendre(order)
|
||||
|
||||
# compute initial guess from previous rho
|
||||
(u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
|
||||
for j in 2:nlrho_init
|
||||
rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
|
||||
(u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
|
||||
end
|
||||
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
@printf("% .15e % .15e\n",x,real(easyeq_u_x(x,2*u-rho*u.*u,weights)))
|
||||
end
|
||||
end
|
||||
|
||||
# condensate fraction
|
||||
function easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx)
|
||||
# compute gaussian quadrature weights
|
||||
weights=gausslegendre(order)
|
||||
|
||||
# compute initial guess from previous rho
|
||||
(u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx)
|
||||
for j in 2:nlrho_init
|
||||
rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1))
|
||||
(u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx)
|
||||
end
|
||||
|
||||
# compute eta
|
||||
eta=easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx)
|
||||
|
||||
# print energy
|
||||
@printf("% .15e % .15e\n",eta,err)
|
||||
end
|
||||
|
||||
# condensate fraction as a function of rho
|
||||
function easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,approx)
|
||||
weights=gausslegendre(order)
|
||||
# init u
|
||||
u=easyeq_init_u(a0,order,weights)
|
||||
|
||||
for j in 1:length(rhos)
|
||||
# compute u (init newton with previously computed u)
|
||||
(u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx)
|
||||
|
||||
# compute eta
|
||||
eta=easyeq_eta(u,order,rhos[j],v,maxiter,tolerance,weights,approx)
|
||||
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],eta,err)
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# initialize u
|
||||
@everywhere function easyeq_init_u(a0,order,weights)
|
||||
u=zeros(Float64,order)
|
||||
for j in 1:order
|
||||
# transformed k
|
||||
k=(1-weights[1][j])/(1+weights[1][j])
|
||||
u[j]=4*pi*a0/k^2
|
||||
end
|
||||
|
||||
return u
|
||||
end
|
||||
|
||||
# \hat u(k) computed using Newton
|
||||
@everywhere function easyeq_hatu(u0,order,rho,v,maxiter,tolerance,weights,approx)
|
||||
# initialize V and Eta
|
||||
(V,V0)=easyeq_init_v(weights,v)
|
||||
(Eta,Eta0)=easyeq_init_H(weights,v)
|
||||
|
||||
# init u
|
||||
u=rho*u0
|
||||
|
||||
# iterate
|
||||
err=Inf
|
||||
for i in 1:maxiter-1
|
||||
new=u-inv(easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx)
|
||||
|
||||
err=norm(new-u)/norm(u)
|
||||
if(err<tolerance)
|
||||
u=new
|
||||
break
|
||||
end
|
||||
u=new
|
||||
end
|
||||
|
||||
return (u/rho,easyeq_en(u,V0,Eta0,rho,weights)*rho/2,err)
|
||||
end
|
||||
|
||||
# \Eta
|
||||
@everywhere function easyeq_H(x,t,weights,v)
|
||||
return (x>t ? 2*t/x : 2)* integrate_legendre(y->2*pi*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y)),0,1,weights)
|
||||
end
|
||||
|
||||
# initialize V
|
||||
@everywhere function easyeq_init_v(weights,v)
|
||||
order=length(weights[1])
|
||||
V=Array{Float64}(undef,order)
|
||||
V0=v(0)
|
||||
for i in 1:order
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
V[i]=v(k)
|
||||
end
|
||||
return(V,V0)
|
||||
end
|
||||
|
||||
# initialize Eta
|
||||
@everywhere function easyeq_init_H(weights,v)
|
||||
order=length(weights[1])
|
||||
Eta=Array{Array{Float64}}(undef,order)
|
||||
Eta0=Array{Float64}(undef,order)
|
||||
for i in 1:order
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
Eta[i]=Array{Float64}(undef,order)
|
||||
for j in 1:order
|
||||
y=(weights[1][j]+1)/2
|
||||
Eta[i][j]=easyeq_H(k,(1-y)/y,weights,v)
|
||||
end
|
||||
y=(weights[1][i]+1)/2
|
||||
Eta0[i]=easyeq_H(0,(1-y)/y,weights,v)
|
||||
end
|
||||
return(Eta,Eta0)
|
||||
end
|
||||
|
||||
# Xi(u)
|
||||
@everywhere function easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx)
|
||||
order=length(weights[1])
|
||||
|
||||
# init
|
||||
out=zeros(Float64,order)
|
||||
|
||||
# compute E before running the loop
|
||||
E=easyeq_en(u,V0,Eta0,rho,weights)
|
||||
|
||||
for i in 1:order
|
||||
# k_i
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
# S_i
|
||||
S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
|
||||
|
||||
# A_K,i
|
||||
A=0.
|
||||
if approx.bK!=0.
|
||||
A+=approx.bK*S
|
||||
end
|
||||
if approx.bK!=1.
|
||||
A+=(1-approx.bK)*E
|
||||
end
|
||||
|
||||
# T
|
||||
if approx.bK==1.
|
||||
T=1.
|
||||
else
|
||||
T=S/A
|
||||
end
|
||||
|
||||
# B
|
||||
if approx.bK==approx.bL
|
||||
B=1.
|
||||
else
|
||||
B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E))
|
||||
end
|
||||
|
||||
# X_i
|
||||
X=k^2/(2*A*rho)
|
||||
|
||||
# U_i
|
||||
out[i]=u[i]-T/(2*(X+1))*Phi(B*T/(X+1)^2)
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# derivative of Xi
|
||||
@everywhere function easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx)
|
||||
order=length(weights[1])
|
||||
|
||||
# init
|
||||
out=zeros(Float64,order,order)
|
||||
|
||||
# compute E before the loop
|
||||
E=easyeq_en(u,V0,Eta0,rho,weights)
|
||||
|
||||
for i in 1:order
|
||||
# k_i
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
# S_i
|
||||
S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
|
||||
|
||||
# A_K,i
|
||||
A=0.
|
||||
if approx.bK!=0.
|
||||
A+=approx.bK*S
|
||||
end
|
||||
if approx.bK!=1.
|
||||
A+=(1-approx.bK)*E
|
||||
end
|
||||
|
||||
# T
|
||||
if approx.bK==1.
|
||||
T=1.
|
||||
else
|
||||
T=S/A
|
||||
end
|
||||
|
||||
# B
|
||||
if approx.bK==approx.bL
|
||||
B=1.
|
||||
else
|
||||
B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E))
|
||||
end
|
||||
|
||||
# X_i
|
||||
X=k^2/(2*A*rho)
|
||||
|
||||
for j in 1:order
|
||||
y=(weights[1][j]+1)/2
|
||||
dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)*weights[2][j]
|
||||
dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^3*y^3)*weights[2][j]
|
||||
|
||||
# dA
|
||||
dA=0.
|
||||
if approx.bK!=0.
|
||||
dA+=approx.bK*dS
|
||||
end
|
||||
if approx.bK!=1.
|
||||
dA+=(1-approx.bK)*dE
|
||||
end
|
||||
|
||||
# dT
|
||||
if approx.bK==1.
|
||||
dT=0.
|
||||
else
|
||||
dT=(1-approx.bK)*(E*dS-S*dE)/A^2
|
||||
end
|
||||
|
||||
# dB
|
||||
if approx.bK==approx.bL
|
||||
dB=0.
|
||||
else
|
||||
dB=(approx.bL*(1-approx.bK)-approx.bK*(1-approx.bL))*(E*dS-S*dE)/(approx.bK*S+(1-approx.bK*E))^2
|
||||
end
|
||||
|
||||
dX=-k^2/(2*A^2*rho)*dA
|
||||
|
||||
out[i,j]=(i==j ? 1 : 0)-(dT-T*dX/(X+1))/(2*(X+1))*Phi(B*T/(X+1)^2)-T/(2*(X+1)^3)*(B*dT+T*dB-2*B*T*dX/(X+1))*dPhi(B*T/(X+1)^2)
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# derivative of Xi with respect to mu
|
||||
@everywhere function easyeq_dXidmu(u,V,V0,Eta,Eta0,weights,rho,approx)
|
||||
order=length(weights[1])
|
||||
|
||||
# init
|
||||
out=zeros(Float64,order)
|
||||
|
||||
# compute E before running the loop
|
||||
E=easyeq_en(u,V0,Eta0,rho,weights)
|
||||
|
||||
for i in 1:order
|
||||
# k_i
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
# S_i
|
||||
S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights)
|
||||
|
||||
# A_K,i
|
||||
A=0.
|
||||
if approx.bK!=0.
|
||||
A+=approx.bK*S
|
||||
end
|
||||
if approx.bK!=1.
|
||||
A+=(1-approx.bK)*E
|
||||
end
|
||||
|
||||
# T
|
||||
if approx.bK==1.
|
||||
T=1.
|
||||
else
|
||||
T=S/A
|
||||
end
|
||||
|
||||
# B
|
||||
if approx.bK==approx.bL
|
||||
B=1.
|
||||
else
|
||||
B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E))
|
||||
end
|
||||
|
||||
# X_i
|
||||
X=k^2/(2*A*rho)
|
||||
|
||||
out[i]=T/(2*rho*A*(X+1)^2)*Phi(B*T/(X+1)^2)+B*T^2/(rho*A*(X+1)^4)*dPhi(B*T/(X+1)^2)
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# energy
|
||||
@everywhere function easyeq_en(u,V0,Eta0,rho,weights)
|
||||
return V0-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights)
|
||||
end
|
||||
|
||||
# condensate fraction
|
||||
@everywhere function easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx)
|
||||
(V,V0)=easyeq_init_v(weights,v)
|
||||
(Eta,Eta0)=easyeq_init_H(weights,v)
|
||||
|
||||
du=-inv(easyeq_dXi(rho*u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_dXidmu(rho*u,V,V0,Eta,Eta0,weights,rho,approx)
|
||||
|
||||
eta=-1/(2*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0,1,weights)
|
||||
|
||||
return eta
|
||||
end
|
||||
|
||||
# inverse Fourier transform
|
||||
@everywhere function easyeq_u_x(x,u,weights)
|
||||
order=length(weights[1])
|
||||
out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0,1,weights)
|
||||
return out
|
||||
end
|
||||
58
src/integration.jl
Normal file
58
src/integration.jl
Normal file
@ -0,0 +1,58 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# approximate \int_a^b f using Gauss-Legendre quadratures
|
||||
@everywhere function integrate_legendre(f,a,b,weights)
|
||||
out=0
|
||||
for i in 1:length(weights[1])
|
||||
out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)
|
||||
end
|
||||
return out
|
||||
end
|
||||
# \int f*g where g is sampled at the Legendre nodes
|
||||
@everywhere function integrate_legendre_sampled(f,g,a,b,weights)
|
||||
out=0
|
||||
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]
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
|
||||
# approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures
|
||||
@everywhere function integrate_chebyshev(f,a,b,N)
|
||||
out=0
|
||||
for i in 1:N
|
||||
out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2)
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# approximate \int_0^\infty dr f(r)*exp(-a*r) using Gauss-Chebyshev quadratures
|
||||
@everywhere function integrate_laguerre(f,a,weights_gL)
|
||||
out=0.
|
||||
for i in 1:length(weights_gL[1])
|
||||
out+=1/a*f(weights_gL[1][i]/a)*weights_gL[2][i]
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# Hann window
|
||||
@everywhere function hann(x,L)
|
||||
if abs(x)<L/2
|
||||
return cos(pi*x/L)^2
|
||||
else
|
||||
return 0.
|
||||
end
|
||||
end
|
||||
108
src/interpolation.jl
Normal file
108
src/interpolation.jl
Normal file
@ -0,0 +1,108 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# linear interpolation: given vectors x,y, compute a linear interpolation for y(x0)
|
||||
# assume x is ordered
|
||||
@everywhere function linear_interpolation(x0,x,y)
|
||||
# if x0 is beyond all x's, then return the corresponding boundary value.
|
||||
if x0>x[length(x)]
|
||||
return y[length(y)]
|
||||
elseif x0<x[1]
|
||||
return y[1]
|
||||
end
|
||||
|
||||
# find bracketing interval
|
||||
i=bracket(x0,x,1,length(x))
|
||||
|
||||
# interpolate
|
||||
return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i])
|
||||
end
|
||||
@everywhere function bracket(x0,x,a,b)
|
||||
i=floor(Int64,(a+b)/2)
|
||||
if x0<x[i]
|
||||
return bracket(x0,x,a,i)
|
||||
elseif x0>x[i+1]
|
||||
return bracket(x0,x,i+1,b)
|
||||
else
|
||||
return i
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# polynomial interpolation of a family of points
|
||||
@everywhere function poly_interpolation(x,y)
|
||||
# init for recursion
|
||||
rec=Array{Polynomial{Float64}}(undef,length(x))
|
||||
for i in 1:length(x)
|
||||
rec[i]=Polynomial([1.])
|
||||
end
|
||||
|
||||
# compute \prod (x-x_i)
|
||||
poly_interpolation_rec(rec,x,1,length(x))
|
||||
|
||||
# sum terms together
|
||||
out=0.
|
||||
for i in 1:length(y)
|
||||
out+=rec[i]/rec[i](x[i])*y[i]
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
# recursive helper function
|
||||
@everywhere function poly_interpolation_rec(out,x,a,b)
|
||||
if a==b
|
||||
return
|
||||
end
|
||||
# mid point
|
||||
mid=floor(Int64,(a+b)/2)
|
||||
# multiply left and right of mid
|
||||
right=Polynomial([1.])
|
||||
for i in mid+1:b
|
||||
right*=Polynomial([-x[i],1.])
|
||||
end
|
||||
left=Polynomial([1.])
|
||||
for i in a:mid
|
||||
left*=Polynomial([-x[i],1.])
|
||||
end
|
||||
|
||||
# multiply into left and right
|
||||
for i in a:mid
|
||||
out[i]*=right
|
||||
end
|
||||
for i in mid+1:b
|
||||
out[i]*=left
|
||||
end
|
||||
|
||||
# recurse
|
||||
poly_interpolation_rec(out,x,a,mid)
|
||||
poly_interpolation_rec(out,x,mid+1,b)
|
||||
|
||||
return
|
||||
end
|
||||
## the following does the same, but has complexity N^2, the function above has N*log(N)
|
||||
#@everywhere function poly_interpolation(x,y)
|
||||
# out=Polynomial([0.])
|
||||
# for i in 1:length(x)
|
||||
# prod=Polynomial([1.])
|
||||
# for j in 1:length(x)
|
||||
# if j!=i
|
||||
# prod*=Polynomial([-x[j]/(x[i]-x[j]),1/(x[i]-x[j])])
|
||||
# end
|
||||
# end
|
||||
# out+=prod*y[i]
|
||||
# end
|
||||
#
|
||||
# return out
|
||||
#end
|
||||
|
||||
443
src/main.jl
Normal file
443
src/main.jl
Normal file
@ -0,0 +1,443 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
using Distributed
|
||||
@everywhere using FastGaussQuadrature
|
||||
@everywhere using LinearAlgebra
|
||||
@everywhere using Polynomials
|
||||
@everywhere using Printf
|
||||
@everywhere using SpecialFunctions
|
||||
include("chebyshev.jl")
|
||||
include("integration.jl")
|
||||
include("interpolation.jl")
|
||||
include("tools.jl")
|
||||
include("potentials.jl")
|
||||
include("print.jl")
|
||||
include("easyeq.jl")
|
||||
include("simpleq-Kv.jl")
|
||||
include("simpleq-iteration.jl")
|
||||
include("simpleq-hardcore.jl")
|
||||
include("anyeq.jl")
|
||||
function main()
|
||||
## defaults
|
||||
|
||||
# values of rho,e, when needed
|
||||
rho=1e-6
|
||||
e=1e-4
|
||||
|
||||
# incrementally initialize Newton algorithm
|
||||
nlrho_init=1
|
||||
|
||||
# potential
|
||||
v=k->v_exp(k,1.)
|
||||
a0=a0_exp(1.)
|
||||
# arguments of the potential
|
||||
v_param_a=1.
|
||||
v_param_b=1.
|
||||
v_param_c=1.
|
||||
v_param_e=1.
|
||||
|
||||
# plot range when plotting in rho
|
||||
minlrho=-6
|
||||
maxlrho=2
|
||||
nlrho=100
|
||||
rhos=Array{Float64}(undef,0)
|
||||
# plot range when plotting in e
|
||||
minle=-6
|
||||
maxle=2
|
||||
nle=100
|
||||
es=Array{Float64}(undef,0)
|
||||
# plot range when plotting in x
|
||||
xmin=0
|
||||
xmax=100
|
||||
nx=100
|
||||
|
||||
# cutoffs
|
||||
tolerance=1e-11
|
||||
order=100
|
||||
maxiter=21
|
||||
|
||||
# for anyeq
|
||||
# P
|
||||
P=11
|
||||
# N
|
||||
N=12
|
||||
# number of splines
|
||||
J=10
|
||||
|
||||
# starting rho from which to incrementally initialize Newton algorithm
|
||||
# default must be set after reading rho, if not set explicitly
|
||||
minlrho_init=nothing
|
||||
|
||||
# Hann window for Fourier transforms
|
||||
windowL=1e3
|
||||
|
||||
# approximation for easyeq
|
||||
# bK,bL
|
||||
easyeq_simpleq_approx=Easyeq_approx(0.,0.)
|
||||
easyeq_medeq_approx=Easyeq_approx(1.,1.)
|
||||
easyeq_approx=easyeq_simpleq_approx
|
||||
|
||||
# approximation for anyeq
|
||||
# aK,bK,gK,aL1,bL1,aL2,bL2,gL2,aL3,bL3,gl3
|
||||
anyeq_simpleq_approx=Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)
|
||||
anyeq_bigeq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,0.,0.,0.)
|
||||
anyeq_fulleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,0.,1.)
|
||||
anyeq_medeq_approx=Anyeq_approx(0.,1.,1.,0.,1.,0.,0.,0.,0.,0.,0.)
|
||||
anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.)
|
||||
anyeq_approx=anyeq_bigeq_approx
|
||||
|
||||
# read cli arguments
|
||||
(params,potential,method,savefile,command)=read_args(ARGS)
|
||||
|
||||
# read params
|
||||
if params!=""
|
||||
for param in split(params,";")
|
||||
terms=split(param,"=")
|
||||
if length(terms)!=2
|
||||
print(stderr,"error: could not read parameter '",param,"'.\n")
|
||||
exit(-1)
|
||||
end
|
||||
lhs=terms[1]
|
||||
rhs=terms[2]
|
||||
if lhs=="rho"
|
||||
rho=parse(Float64,rhs)
|
||||
elseif lhs=="minlrho_init"
|
||||
minlrho_init=parse(Float64,rhs)
|
||||
elseif lhs=="nlrho_init"
|
||||
nlrho_init=parse(Int64,rhs)
|
||||
elseif lhs=="e"
|
||||
e=parse(Float64,rhs)
|
||||
elseif lhs=="tolerance"
|
||||
tolerance=parse(Float64,rhs)
|
||||
elseif lhs=="order"
|
||||
order=parse(Int64,rhs)
|
||||
elseif lhs=="maxiter"
|
||||
maxiter=parse(Int64,rhs)
|
||||
elseif lhs=="rhos"
|
||||
rhos=parse_list(rhs)
|
||||
elseif lhs=="minlrho"
|
||||
minlrho=parse(Float64,rhs)
|
||||
elseif lhs=="maxlrho"
|
||||
maxlrho=parse(Float64,rhs)
|
||||
elseif lhs=="nlrho"
|
||||
nlrho=parse(Int64,rhs)
|
||||
elseif lhs=="es"
|
||||
es=parse_list(rhs)
|
||||
elseif lhs=="minle"
|
||||
minle=parse(Float64,rhs)
|
||||
elseif lhs=="maxle"
|
||||
maxle=parse(Float64,rhs)
|
||||
elseif lhs=="nle"
|
||||
nle=parse(Int64,rhs)
|
||||
elseif lhs=="xmin"
|
||||
xmin=parse(Float64,rhs)
|
||||
elseif lhs=="xmax"
|
||||
xmax=parse(Float64,rhs)
|
||||
elseif lhs=="nx"
|
||||
nx=parse(Int64,rhs)
|
||||
elseif lhs=="P"
|
||||
P=parse(Int64,rhs)
|
||||
elseif lhs=="N"
|
||||
N=parse(Int64,rhs)
|
||||
elseif lhs=="J"
|
||||
J=parse(Int64,rhs)
|
||||
elseif lhs=="window_L"
|
||||
windowL=parse(Float64,rhs)
|
||||
elseif lhs=="aK"
|
||||
anyeq_approx.aK=parse(Float64,rhs)
|
||||
elseif lhs=="bK"
|
||||
anyeq_approx.bK=parse(Float64,rhs)
|
||||
easyeq_approx.bK=parse(Float64,rhs)
|
||||
elseif lhs=="gK"
|
||||
anyeq_approx.gK=parse(Float64,rhs)
|
||||
elseif lhs=="aL1"
|
||||
anyeq_approx.aL1=parse(Float64,rhs)
|
||||
elseif lhs=="bL"
|
||||
easyeq_approx.bL=parse(Float64,rhs)
|
||||
elseif lhs=="bL1"
|
||||
anyeq_approx.bL1=parse(Float64,rhs)
|
||||
elseif lhs=="aL2"
|
||||
anyeq_approx.aL2=parse(Float64,rhs)
|
||||
elseif lhs=="bL2"
|
||||
anyeq_approx.bL2=parse(Float64,rhs)
|
||||
elseif lhs=="gL2"
|
||||
anyeq_approx.gL2=parse(Float64,rhs)
|
||||
elseif lhs=="aL3"
|
||||
anyeq_approx.aL3=parse(Float64,rhs)
|
||||
elseif lhs=="bL3"
|
||||
anyeq_approx.bL3=parse(Float64,rhs)
|
||||
elseif lhs=="gL3"
|
||||
anyeq_approx.gL3=parse(Float64,rhs)
|
||||
elseif lhs=="v_a"
|
||||
v_param_a=parse(Float64,rhs)
|
||||
elseif lhs=="v_b"
|
||||
v_param_b=parse(Float64,rhs)
|
||||
elseif lhs=="v_c"
|
||||
v_param_c=parse(Float64,rhs)
|
||||
elseif lhs=="v_e"
|
||||
v_param_e=parse(Float64,rhs)
|
||||
elseif lhs=="eq"
|
||||
if rhs=="simpleq"
|
||||
easyeq_approx=easyeq_simpleq_approx
|
||||
anyeq_approx=anyeq_simpleq_approx
|
||||
elseif rhs=="medeq"
|
||||
easyeq_approx=easyeq_medeq_approx
|
||||
anyeq_approx=anyeq_medeq_approx
|
||||
elseif rhs=="bigeq"
|
||||
anyeq_approx=anyeq_bigeq_approx
|
||||
elseif rhs=="fulleq"
|
||||
anyeq_approx=anyeq_fulleq_approx
|
||||
elseif rhs=="compleq"
|
||||
anyeq_approx=anyeq_compleq_approx
|
||||
else
|
||||
print(stderr,"error: unrecognized equation: ",rhs,"\n")
|
||||
exit(-1)
|
||||
end
|
||||
else
|
||||
print(stderr,"unrecognized parameter '",lhs,"'.\n")
|
||||
exit(-1)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
## read potential
|
||||
if potential=="exp"
|
||||
v=k->v_exp(k,v_param_a)
|
||||
a0=a0_exp(v_param_a)
|
||||
elseif potential=="expcry"
|
||||
v=k->v_expcry(k,v_param_a,v_param_b)
|
||||
a0=a0_expcry(v_param_a,v_param_b)
|
||||
elseif potential=="npt"
|
||||
v=k->v_npt(k)
|
||||
a0=a0_npt()
|
||||
elseif potential=="alg"
|
||||
v=v_alg
|
||||
a0=a0_alg
|
||||
elseif potential=="algwell"
|
||||
v=v_algwell
|
||||
a0=a0_algwell
|
||||
elseif potential=="exact"
|
||||
v=k->v_exact(k,v_param_b,v_param_c,v_param_e)
|
||||
a0=a0_exact(v_param_b,v_param_c,v_param_e)
|
||||
elseif potential=="tent"
|
||||
v=k->v_tent(k,v_param_a,v_param_b)
|
||||
a0=a0_tent(v_param_a,v_param_b)
|
||||
else
|
||||
print(stderr,"unrecognized potential '",potential,"'.\n'")
|
||||
exit(-1)
|
||||
end
|
||||
|
||||
## set parameters
|
||||
# rhos
|
||||
if length(rhos)==0
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
end
|
||||
# es
|
||||
if length(es)==0
|
||||
es=Array{Float64}(undef,nle)
|
||||
for j in 0:nle-1
|
||||
es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j))
|
||||
end
|
||||
end
|
||||
|
||||
# default minlrho_init
|
||||
if (minlrho_init==nothing)
|
||||
minlrho_init=log10(rho)
|
||||
end
|
||||
|
||||
# splines
|
||||
taus=Array{Float64}(undef,J+1)
|
||||
for j in 0:J
|
||||
taus[j+1]=-1+2*j/J
|
||||
end
|
||||
|
||||
## run command
|
||||
if method=="easyeq"
|
||||
if command=="energy"
|
||||
easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
|
||||
# e(rho)
|
||||
elseif command=="energy_rho"
|
||||
easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx)
|
||||
# u(k)
|
||||
elseif command=="uk"
|
||||
easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
|
||||
# u(x)
|
||||
elseif command=="ux"
|
||||
easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,easyeq_approx)
|
||||
elseif command=="uux"
|
||||
easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,easyeq_approx)
|
||||
# condensate fraction
|
||||
elseif command=="condensate_fraction"
|
||||
easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx)
|
||||
elseif command=="condensate_fraction_rho"
|
||||
easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx)
|
||||
else
|
||||
print(stderr,"unrecognized command '",command,"'.\n")
|
||||
exit(-1)
|
||||
end
|
||||
elseif method=="simpleq-hardcore"
|
||||
if command=="energy_rho"
|
||||
simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
|
||||
elseif command=="ux"
|
||||
simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
|
||||
elseif command=="condensate_fraction_rho"
|
||||
simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance)
|
||||
end
|
||||
elseif method=="simpleq-iteration"
|
||||
# u_n(x) using iteration
|
||||
if command=="u"
|
||||
simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx)
|
||||
# rho(e) using iteration
|
||||
elseif command=="rho_e"
|
||||
simpleq_iteration_rho_e(es,order,v,maxiter)
|
||||
else
|
||||
print(stderr,"unrecognized command '",command,"'.\n")
|
||||
exit(-1)
|
||||
end
|
||||
elseif method=="anyeq"
|
||||
if command=="save_Abar"
|
||||
anyeq_save_Abar(taus,P,N,J,v,anyeq_approx)
|
||||
elseif command=="energy"
|
||||
anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
|
||||
# e(rho)
|
||||
elseif command=="energy_rho"
|
||||
anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
|
||||
elseif command=="energy_rho_init_prevrho"
|
||||
anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
|
||||
elseif command=="energy_rho_init_nextrho"
|
||||
anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
|
||||
# u(j)
|
||||
elseif command=="uk"
|
||||
anyeq_uk(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,anyeq_approx,savefile)
|
||||
# u(j)
|
||||
elseif command=="ux"
|
||||
anyeq_ux(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile)
|
||||
# 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)
|
||||
elseif command=="condensate_fraction_rho"
|
||||
anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile)
|
||||
# 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)
|
||||
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)
|
||||
else
|
||||
print(stderr,"unrecognized command: '",command,"'\n")
|
||||
exit(-1)
|
||||
end
|
||||
elseif method=="simpleq-Kv"
|
||||
if command=="2pt"
|
||||
simpleq_Kv_2pt(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx)
|
||||
elseif command=="condensate_fraction_rho"
|
||||
simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
|
||||
elseif command=="Kv"
|
||||
simpleq_Kv_Kv(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx)
|
||||
else
|
||||
print(stderr,"unrecognized command: '",command,"'\n")
|
||||
exit(-1)
|
||||
end
|
||||
else
|
||||
print(stderr,"unrecognized method: '",method,"'\n")
|
||||
exit(-1)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
# parse a comma separated list as an array of Float64
|
||||
function parse_list(str)
|
||||
elems=split(str,",")
|
||||
out=Array{Float64}(undef,length(elems))
|
||||
for i in 1:length(elems)
|
||||
out[i]=parse(Float64,elems[i])
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# read cli arguments
|
||||
function read_args(ARGS)
|
||||
# flag
|
||||
flag=""
|
||||
|
||||
# output strings
|
||||
params=""
|
||||
# default potential
|
||||
potential="exp"
|
||||
# default method
|
||||
method="easyeq"
|
||||
|
||||
savefile=""
|
||||
command=""
|
||||
|
||||
# loop over arguments
|
||||
for arg in ARGS
|
||||
# argument that starts with a dash
|
||||
if arg[1]=='-'
|
||||
# go through characters after dash
|
||||
for char in arg[2:length(arg)]
|
||||
|
||||
# set params
|
||||
if char=='p'
|
||||
# raise flag
|
||||
flag="params"
|
||||
elseif char=='U'
|
||||
# raise flag
|
||||
flag="potential"
|
||||
elseif char=='M'
|
||||
# raise flag
|
||||
flag="method"
|
||||
elseif char=='s'
|
||||
# raise flag
|
||||
flag="savefile"
|
||||
else
|
||||
print_usage()
|
||||
exit(-1)
|
||||
end
|
||||
end
|
||||
# arguments that do not start with a dash
|
||||
else
|
||||
if flag=="params"
|
||||
params=arg
|
||||
elseif flag=="potential"
|
||||
potential=arg
|
||||
elseif flag=="method"
|
||||
method=arg
|
||||
elseif flag=="savefile"
|
||||
savefile=arg
|
||||
else
|
||||
command=arg
|
||||
end
|
||||
# reset flag
|
||||
flag=""
|
||||
end
|
||||
end
|
||||
|
||||
if command==""
|
||||
print_usage()
|
||||
exit(-1)
|
||||
end
|
||||
|
||||
return (params,potential,method,savefile,command)
|
||||
end
|
||||
|
||||
# usage
|
||||
function print_usage()
|
||||
print(stderr,"usage: simplesolv [-p params] [-U potential] [-M method] [-s savefile] <command>\n")
|
||||
end
|
||||
|
||||
main()
|
||||
84
src/potentials.jl
Normal file
84
src/potentials.jl
Normal file
@ -0,0 +1,84 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# exponential potential in 3 dimensions
|
||||
@everywhere function v_exp(k,a)
|
||||
return 8*pi/(1+k^2)^2*a
|
||||
end
|
||||
@everywhere function a0_exp(a)
|
||||
if a>0.
|
||||
return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a))
|
||||
elseif a<0.
|
||||
return log(-a)+2*MathConstants.eulergamma-pi*bessely(0,2*sqrt(-a))/besselj(0,2*sqrt(-a))
|
||||
else
|
||||
return 0.
|
||||
end
|
||||
end
|
||||
|
||||
# exp(-x)-a*exp(-b*x) in 3 dimensions
|
||||
@everywhere function v_expcry(k,a,b)
|
||||
return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2))
|
||||
end
|
||||
@everywhere function a0_expcry(a,b)
|
||||
return 1.21751642717932720441274114683413710125487579284827462 #ish
|
||||
end
|
||||
|
||||
# x^2*exp(-|x|) in 3 dimensions
|
||||
@everywhere function v_npt(k)
|
||||
return 96*pi*(1-k^2)/(1+k^2)^4
|
||||
end
|
||||
@everywhere function a0_npt()
|
||||
return 1. #ish
|
||||
end
|
||||
|
||||
# 1/(1+x^4/4) potential in 3 dimensions
|
||||
@everywhere function v_alg(k)
|
||||
if(k==0)
|
||||
return 4*pi^2
|
||||
else
|
||||
return 4*pi^2*exp(-k)*sin(k)/k
|
||||
end
|
||||
end
|
||||
a0_alg=1. #ish
|
||||
|
||||
# (1+a x^4)/(1+x^2)^4 potential in 3 dimensions
|
||||
@everywhere function v_algwell(k)
|
||||
a=4
|
||||
return pi^2/24*exp(-k)*(a*(k^2-9*k+15)+k^2+3*k+3)
|
||||
end
|
||||
a0_algwell=1. #ish
|
||||
|
||||
# potential corresponding to the exact solution c/(1+b^2x^2)^2
|
||||
@everywhere function v_exact(k,b,c,e)
|
||||
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
|
||||
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))
|
||||
end
|
||||
end
|
||||
@everywhere function a0_exact(b,c,e)
|
||||
return 1. #ish
|
||||
end
|
||||
|
||||
# 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)
|
||||
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
|
||||
else
|
||||
return (b/2)^3*a*(4*pi/3)^2
|
||||
end
|
||||
end
|
||||
@everywhere function a0_tent(a,b)
|
||||
return b #ish
|
||||
end
|
||||
52
src/print.jl
Normal file
52
src/print.jl
Normal file
@ -0,0 +1,52 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# print progress
|
||||
@everywhere function progress(
|
||||
j::Int64,
|
||||
tot::Int64,
|
||||
freq::Int64
|
||||
)
|
||||
if (j-1)%ceil(Int,tot/freq)==0
|
||||
if j>1
|
||||
@printf(stderr,"\r")
|
||||
end
|
||||
@printf(stderr,"%d/%d",j,tot)
|
||||
end
|
||||
if j==tot
|
||||
@printf(stderr,"\r")
|
||||
@printf(stderr,"%d/%d\n",j,tot)
|
||||
end
|
||||
end
|
||||
|
||||
# print progress of two indices at once
|
||||
@everywhere function progress_mat(
|
||||
j1::Int64,
|
||||
tot1::Int64,
|
||||
j2::Int64,
|
||||
tot2::Int64,
|
||||
freq::Int64
|
||||
)
|
||||
if ((j1-1)*tot2+j2-1)%ceil(Int,tot1*tot2/freq)==0
|
||||
if j1>1 || j2>1
|
||||
@printf(stderr,"\r")
|
||||
end
|
||||
@printf(stderr,"%2d/%2d, %2d/%2d",j1,tot1,j2,tot2)
|
||||
end
|
||||
if j1==tot1 && j2==tot2
|
||||
@printf(stderr,"\r")
|
||||
@printf(stderr,"%2d/%2d, %2d/%2d\n",j1,tot1,j2,tot2)
|
||||
end
|
||||
end
|
||||
|
||||
119
src/simpleq-Kv.jl
Normal file
119
src/simpleq-Kv.jl
Normal file
@ -0,0 +1,119 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# 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)
|
||||
# init vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
u0=u0s[nlrho]
|
||||
|
||||
(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
|
||||
Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
|
||||
|
||||
# switch to real space
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
Kv=0.
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
Kv+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*Kvk[zeta*N+j]*sin(k[zeta*N+j]*x)
|
||||
end
|
||||
end
|
||||
@printf("% .15e % .15e\n",x,Kv)
|
||||
end
|
||||
end
|
||||
|
||||
# Compute the condensate fraction for simpleq using Kv
|
||||
function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance)
|
||||
# init vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# compute initial guess from medeq
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,0.,maxiter,tolerance)
|
||||
|
||||
for j in 1:length(rhos)
|
||||
(u,E,error)=anyeq_hatu(u0s[j],0.,P,N,J,rhos[j],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
|
||||
Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
|
||||
|
||||
# denominator: (1-rho*\int (Kv)(2u-rho u*u))
|
||||
denom=1-rhos[j]*integrate_f_chebyshev(s->1.,Kvk.*(2*u-rhos[j]*u.^2),k,taus,weights,N,J)
|
||||
|
||||
eta=rhos[j]*integrate_f_chebyshev(s->1.,Kvk.*u,k,taus,weights,N,J)/denom
|
||||
|
||||
@printf("% .15e % .15e\n",rhos[j],eta)
|
||||
end
|
||||
end
|
||||
|
||||
# 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)
|
||||
# init vectors
|
||||
(weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v)
|
||||
|
||||
# compute initial guess from medeq
|
||||
rhos=Array{Float64}(undef,nlrho)
|
||||
for j in 0:nlrho-1
|
||||
rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j))
|
||||
end
|
||||
u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance)
|
||||
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.))
|
||||
|
||||
# Kv in Fourier space
|
||||
Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
|
||||
|
||||
# denominator: (1-rho*\int (Kv)(2u-rho u*u))
|
||||
denom=1-rho*integrate_f_chebyshev(s->1.,Kvk.*(2*u-rho*u.^2),k,taus,weights,N,J)
|
||||
|
||||
# Kv, u, u*Kv, u*u*Kv in real space
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
ux=inverse_fourier_chebyshev(u,x,k,taus,weights,N,J)
|
||||
Kv=inverse_fourier_chebyshev(Kvk,x,k,taus,weights,N,J)
|
||||
uKv=inverse_fourier_chebyshev(u.*Kvk,x,k,taus,weights,N,J)
|
||||
uuKv=inverse_fourier_chebyshev(u.*u.*Kvk,x,k,taus,weights,N,J)
|
||||
|
||||
# correlation in real space
|
||||
Cx=rho^2*((1-ux)-((1-ux)*Kv-2*rho*uKv+rho^2*uuKv)/denom)
|
||||
|
||||
@printf("% .15e % .15e\n",x,Cx)
|
||||
end
|
||||
end
|
||||
|
||||
# Kv
|
||||
function simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J)
|
||||
# (-Delta+v+4e(1-\rho u*)) in Fourier space
|
||||
M=Array{Float64,2}(undef,N*J,N*J)
|
||||
for zetapp in 0:J-1
|
||||
for n in 1:N
|
||||
M[:,zetapp*N+n]=conv_one_v_chebyshev(n,zetapp,Upsilon,k,taus,weights,N,J)
|
||||
end
|
||||
end
|
||||
for i in 1:J*N
|
||||
M[i,i]+=k[i]^2+4*E*(1-rho*u[i])
|
||||
end
|
||||
|
||||
return inv(M)*V
|
||||
end
|
||||
573
src/simpleq-hardcore.jl
Normal file
573
src/simpleq-hardcore.jl
Normal file
@ -0,0 +1,573 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# compute energy as a function of rho
|
||||
function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance)
|
||||
## spawn workers
|
||||
# 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]=zeros(0)
|
||||
end
|
||||
for j in 1:length(rhos)
|
||||
append!(work[j%nw+1],j)
|
||||
end
|
||||
|
||||
# initialize vectors
|
||||
(weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
|
||||
|
||||
# initial guess
|
||||
u0s=Array{Array{Float64}}(undef,length(rhos))
|
||||
e0s=Array{Float64}(undef,length(rhos))
|
||||
for j in 1:length(rhos)
|
||||
u0s[j]=Array{Float64}(undef,N*J)
|
||||
for i in 1:N*J
|
||||
u0s[j][i]=1/(1+r[i]^2)^2
|
||||
end
|
||||
e0s[j]=pi*rhos[j]
|
||||
end
|
||||
|
||||
|
||||
# save result from each task
|
||||
us=Array{Array{Float64}}(undef,length(rhos))
|
||||
es=Array{Float64}(undef,length(rhos))
|
||||
err=Array{Float64}(undef,length(rhos))
|
||||
|
||||
count=0
|
||||
# for each worker
|
||||
@sync for p in 1:nw
|
||||
# for each task
|
||||
@async for j in work[p]
|
||||
count=count+1
|
||||
if count>=nw
|
||||
progress(count,length(rhos),10000)
|
||||
end
|
||||
# run the task
|
||||
(us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
|
||||
end
|
||||
end
|
||||
|
||||
for j in 1:length(rhos)
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j])
|
||||
end
|
||||
end
|
||||
|
||||
# compute u(x)
|
||||
function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance)
|
||||
# initialize vectors
|
||||
(weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
|
||||
|
||||
# initial guess
|
||||
u0=Array{Float64}(undef,N*J)
|
||||
for i in 1:N*J
|
||||
u0[i]=1/(1+r[i]^2)^2
|
||||
end
|
||||
e0=pi*rho
|
||||
|
||||
(u,e,err)=simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
|
||||
|
||||
for i in 1:N*J
|
||||
@printf("% .15e % .15e\n",r[i],u[i])
|
||||
end
|
||||
end
|
||||
|
||||
# compute condensate fraction as a function of rho
|
||||
function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance)
|
||||
## spawn workers
|
||||
# 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]=zeros(0)
|
||||
end
|
||||
for j in 1:length(rhos)
|
||||
append!(work[j%nw+1],j)
|
||||
end
|
||||
|
||||
# initialize vectors
|
||||
(weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J)
|
||||
|
||||
# initial guess
|
||||
u0s=Array{Array{Float64}}(undef,length(rhos))
|
||||
e0s=Array{Float64}(undef,length(rhos))
|
||||
for j in 1:length(rhos)
|
||||
u0s[j]=Array{Float64}(undef,N*J)
|
||||
for i in 1:N*J
|
||||
u0s[j][i]=1/(1+r[i]^2)^2
|
||||
end
|
||||
e0s[j]=pi*rhos[j]
|
||||
end
|
||||
|
||||
|
||||
# save result from each task
|
||||
us=Array{Array{Float64}}(undef,length(rhos))
|
||||
es=Array{Float64}(undef,length(rhos))
|
||||
err=Array{Float64}(undef,length(rhos))
|
||||
|
||||
count=0
|
||||
# for each worker
|
||||
@sync for p in 1:nw
|
||||
# for each task
|
||||
@async for j in work[p]
|
||||
count=count+1
|
||||
if count>=nw
|
||||
progress(count,length(rhos),10000)
|
||||
end
|
||||
# run the task
|
||||
(us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
|
||||
end
|
||||
end
|
||||
|
||||
for j in 1:length(rhos)
|
||||
du=-inv(simpleq_hardcore_DXi(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_dXidmu(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J)
|
||||
@printf("% .15e % .15e % .15e\n",rhos[j],du[N*J+1],err[j])
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# initialize computation
|
||||
@everywhere function simpleq_hardcore_init(taus,P,N,J)
|
||||
# Gauss-Legendre weights
|
||||
weights=gausslegendre(N)
|
||||
weights_gL=gausslaguerre(N)
|
||||
|
||||
# r
|
||||
r=Array{Float64}(undef,J*N)
|
||||
for zeta in 0:J-1
|
||||
for j in 1:N
|
||||
xj=weights[1][j]
|
||||
# set kj
|
||||
r[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1])
|
||||
end
|
||||
end
|
||||
|
||||
# Chebyshev polynomials
|
||||
T=chebyshev_polynomials(P)
|
||||
|
||||
return (weights,weights_gL,r,T)
|
||||
end
|
||||
|
||||
# compute u using chebyshev expansions
|
||||
@everywhere function simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance)
|
||||
# init
|
||||
vec=Array{Float64}(undef,J*N+1)
|
||||
for i in 1:J*N
|
||||
vec[i]=u0[i]
|
||||
end
|
||||
vec[J*N+1]=e0
|
||||
|
||||
# quantify relative error
|
||||
error=Inf
|
||||
|
||||
# run Newton algorithm
|
||||
for i in 1:maxiter-1
|
||||
u=vec[1:J*N]
|
||||
e=vec[J*N+1]
|
||||
new=vec-inv(simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
|
||||
|
||||
error=norm(new-vec)/norm(new)
|
||||
if(error<tolerance)
|
||||
vec=new
|
||||
break
|
||||
end
|
||||
|
||||
vec=new
|
||||
end
|
||||
|
||||
return(vec[1:J*N],vec[J*N+1],error)
|
||||
end
|
||||
|
||||
# Xi
|
||||
@everywhere function simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
|
||||
out=Array{Float64}(undef,J*N+1)
|
||||
FU=chebyshev(u,taus,weights,P,N,J,4)
|
||||
|
||||
#D's
|
||||
d0=D0(e,rho,r,weights,N,J)
|
||||
d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
|
||||
# u
|
||||
for i in 1:J*N
|
||||
out[i]=d0[i]+dot(FU,d1[i])+dot(FU,d2[i]*FU)-u[i]
|
||||
end
|
||||
# e
|
||||
out[J*N+1]=-e+2*pi*rho*
|
||||
((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/
|
||||
(1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU)))
|
||||
|
||||
return out
|
||||
end
|
||||
# DXi
|
||||
@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
|
||||
out=Array{Float64,2}(undef,J*N+1,J*N+1)
|
||||
FU=chebyshev(u,taus,weights,P,N,J,4)
|
||||
|
||||
#D's
|
||||
d0=D0(e,rho,r,weights,N,J)
|
||||
d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
dsed0=dseD0(e,rho,r,weights,N,J)
|
||||
dsed1=dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
dsed2=dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
|
||||
# denominator of e
|
||||
denom=1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU))
|
||||
|
||||
for zetapp in 0:J-1
|
||||
for n in 1:N
|
||||
one=zeros(Int64,J*N)
|
||||
one[zetapp*N+n]=1
|
||||
Fone=chebyshev(one,taus,weights,P,N,J,4)
|
||||
|
||||
for i in 1:J*N
|
||||
# du/du
|
||||
out[i,zetapp*N+n]=dot(Fone,d1[i])+2*dot(FU,d2[i]*Fone)-(zetapp*N+n==i ? 1 : 0)
|
||||
# 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)
|
||||
end
|
||||
# de/du
|
||||
out[J*N+1,zetapp*N+n]=2*pi*rho*
|
||||
(-dot(Fone,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-2*dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*Fone))/denom-
|
||||
2*pi*rho*
|
||||
((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))*
|
||||
rho^2*(dot(Fone,gammabar1(taus,T,weights,P,J,4))+2*dot(FU,gammabar2(taus,T,weights,P,J,4)*Fone))/denom^2
|
||||
end
|
||||
end
|
||||
#de/de
|
||||
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*sqrt(abs(e)))*(e>=0 ? 1 : -1)
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# dXi/dmu
|
||||
@everywhere function simpleq_hardcore_dXidmu(u,e,rho,r,taus,T,weights,weights_gL,P,N,J)
|
||||
out=Array{Float64}(undef,J*N+1)
|
||||
FU=chebyshev(u,taus,weights,P,N,J,4)
|
||||
|
||||
#D's
|
||||
dsmud0=dsmuD0(e,rho,r,weights,N,J)
|
||||
dsmud1=dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
dsmud2=dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4)
|
||||
|
||||
# u
|
||||
for i in 1:J*N
|
||||
out[i]=(dsmud0[i]+dot(FU,dsmud1[i])+dot(FU,dsmud2[i]*FU))/(4*sqrt(abs(e)))
|
||||
end
|
||||
# e
|
||||
out[J*N+1]=2*pi*rho*(2/3+1/(2*sqrt(abs(e)))-
|
||||
(dsmudgamma0(e,rho,weights)+dot(FU,dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))+dot(FU,dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/(4*sqrt(abs(e)))
|
||||
)/(1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU)))
|
||||
|
||||
return out
|
||||
end
|
||||
|
||||
# B's
|
||||
@everywhere function B0(r)
|
||||
return pi/12*(r-1)^2*(r+5)
|
||||
end
|
||||
@everywhere function B1(r,zeta,n,taus,T,weights,nu)
|
||||
return (taus[zeta+1]>=(2-r)/r || taus[zeta+2]<=-r/(r+2) ? 0 :
|
||||
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))
|
||||
,min(taus[zeta+2],(2-r)/r),weights)
|
||||
)
|
||||
end
|
||||
@everywhere function B2(r,zeta,n,zetap,m,taus,T,weights,nu)
|
||||
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]))*
|
||||
(taus[zetap+1]>=alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau) || taus[zetap+2]<=alpham(1+r,tau) ? 0 :
|
||||
integrate_legendre(s->
|
||||
1/(1+s)^(3-nu)*T[m+1]((2*s-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1]
|
||||
,alpham(1+r,tau)),min(taus[zetap+2],alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau)),weights)
|
||||
)
|
||||
,taus[zeta+1],taus[zeta+2],weights)
|
||||
end
|
||||
|
||||
# D's
|
||||
@everywhere function D0(e,rho,r,weights,N,J)
|
||||
out=Array{Float64}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
|
||||
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
|
||||
,0,min(1,r[i]),weights)+
|
||||
(r[i]>=1 ? 0 :
|
||||
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)
|
||||
,0,1-r[i],weights)
|
||||
)
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
|
||||
out=Array{Array{Float64}}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=Array{Float64}(undef,(P+1)*J)
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
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
|
||||
,0,r[i],weights)+
|
||||
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)
|
||||
,2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return out
|
||||
end
|
||||
@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
|
||||
out=Array{Array{Float64,2}}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
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->
|
||||
(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)+
|
||||
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)
|
||||
,2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# dD/d sqrt(abs(e))'s
|
||||
@everywhere function dseD0(e,rho,r,weights,N,J)
|
||||
out=Array{Float64}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
|
||||
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)
|
||||
,0,min(1,r[i]),weights)+
|
||||
(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)
|
||||
)
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
|
||||
out=Array{Array{Float64}}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=Array{Float64}(undef,(P+1)*J)
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
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)
|
||||
,0,r[i],weights)+
|
||||
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]))
|
||||
,2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
|
||||
out=Array{Array{Float64,2}}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
for m in 0:P
|
||||
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)
|
||||
,0,r[i],weights)+
|
||||
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]))
|
||||
,2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# dD/d sqrt(abs(e+mu/2))'s
|
||||
@everywhere function dsmuD0(e,rho,r,weights,N,J)
|
||||
out=Array{Float64}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+
|
||||
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)
|
||||
,0,min(1,r[i]),weights)+
|
||||
(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)
|
||||
)
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
|
||||
out=Array{Array{Float64}}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=Array{Float64}(undef,(P+1)*J)
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
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)
|
||||
,0,r[i],weights)+
|
||||
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]))
|
||||
,2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu)
|
||||
out=Array{Array{Float64,2}}(undef,J*N)
|
||||
for i in 1:J*N
|
||||
out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J)
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
for m in 0:P
|
||||
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)
|
||||
,0,r[i],weights)+
|
||||
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]))
|
||||
,2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# gamma's
|
||||
@everywhere function gamma0(e,rho,weights)
|
||||
return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights)
|
||||
end
|
||||
@everywhere function gamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
|
||||
out=Array{Float64}(undef,J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
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)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
|
||||
out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
for m in 0:P
|
||||
out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=2*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# dgamma/d sqrt(abs(e))'s
|
||||
@everywhere function dsedgamma0(e,rho,weights)
|
||||
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
|
||||
@everywhere function dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
|
||||
out=Array{Float64}(undef,J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
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)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
|
||||
out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
for m in 0:P
|
||||
out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# dgamma/d sqrt(e+mu/2)'s
|
||||
@everywhere function dsmudgamma0(e,rho,weights)
|
||||
return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights)
|
||||
end
|
||||
@everywhere function dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu)
|
||||
out=Array{Float64}(undef,J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
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)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu)
|
||||
out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
for m in 0:P
|
||||
out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# \bar gamma's
|
||||
@everywhere function gammabar0(weights)
|
||||
return 4*pi*integrate_legendre(s->s^2*B0(s-1),0,1,weights)
|
||||
end
|
||||
@everywhere function gammabar1(taus,T,weights,P,J,nu)
|
||||
out=Array{Float64}(undef,J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
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)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function gammabar2(taus,T,weights,P,J,nu)
|
||||
out=Array{Float64,2}(undef,J*(P+1),J*(P+1))
|
||||
for zeta in 0:J-1
|
||||
for n in 0:P
|
||||
for zetap in 0:J-1
|
||||
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)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
94
src/simpleq-iteration.jl
Normal file
94
src/simpleq-iteration.jl
Normal file
@ -0,0 +1,94 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# compute rho(e) using the iteration
|
||||
function simpleq_iteration_rho_e(es,order,v,maxiter)
|
||||
for j in 1:length(es)
|
||||
(u,rho)=simpleq_iteration_hatun(es[j],order,v,maxiter)
|
||||
@printf("% .15e % .15e\n",es[j],real(rho[maxiter+1]))
|
||||
end
|
||||
end
|
||||
|
||||
# compute u(x) using the iteration and print at every step
|
||||
function simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx)
|
||||
(u,rho)=simpleq_iteration_hatun(e,order,v,maxiter)
|
||||
|
||||
weights=gausslegendre(order)
|
||||
for i in 1:nx
|
||||
x=xmin+(xmax-xmin)*i/nx
|
||||
@printf("% .15e ",x)
|
||||
for n in 2:maxiter+1
|
||||
@printf("% .15e ",real(easyeq_u_x(x,u[n],weights)))
|
||||
end
|
||||
print('\n')
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
# \hat u_n
|
||||
function simpleq_iteration_hatun(e,order,v,maxiter)
|
||||
# gauss legendre weights
|
||||
weights=gausslegendre(order)
|
||||
|
||||
# initialize V and Eta
|
||||
(V,V0)=easyeq_init_v(weights,v)
|
||||
(Eta,Eta0)=easyeq_init_H(weights,v)
|
||||
|
||||
# init u and rho
|
||||
u=Array{Array{Float64}}(undef,maxiter+1)
|
||||
u[1]=zeros(Float64,order)
|
||||
rho=zeros(Float64,maxiter+1)
|
||||
|
||||
# iterate
|
||||
for n in 1:maxiter
|
||||
u[n+1]=simpleq_iteration_A(e,weights,Eta)\simpleq_iteration_b(u[n],e,rho[n],V)
|
||||
rho[n+1]=simpleq_iteration_rhon(u[n+1],e,weights,V0,Eta0)
|
||||
end
|
||||
|
||||
return (u,rho)
|
||||
end
|
||||
|
||||
# A
|
||||
function simpleq_iteration_A(e,weights,Eta)
|
||||
N=length(weights[1])
|
||||
out=zeros(Float64,N,N)
|
||||
for i in 1:N
|
||||
k=(1-weights[1][i])/(1+weights[1][i])
|
||||
out[i,i]=k^2+4*e
|
||||
for j in 1:N
|
||||
y=(weights[1][j]+1)/2
|
||||
out[i,j]+=weights[2][j]*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)
|
||||
end
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# b
|
||||
function simpleq_iteration_b(u,e,rho,V)
|
||||
out=zeros(Float64,length(V))
|
||||
for i in 1:length(V)
|
||||
out[i]=V[i]+2*e*rho*u[i]^2
|
||||
end
|
||||
return out
|
||||
end
|
||||
|
||||
# rho_n
|
||||
function simpleq_iteration_rhon(u,e,weights,V0,Eta0)
|
||||
S=V0
|
||||
for i in 1:length(weights[1])
|
||||
y=(weights[1][i]+1)/2
|
||||
S+=-weights[2][i]*(1-y)*u[i]*Eta0[i]/(2*(2*pi)^3*y^3)
|
||||
end
|
||||
return 2*e/S
|
||||
end
|
||||
49
src/tools.jl
Normal file
49
src/tools.jl
Normal file
@ -0,0 +1,49 @@
|
||||
## Copyright 2021 Ian Jauslin
|
||||
##
|
||||
## Licensed under the Apache License, Version 2.0 (the "License");
|
||||
## you may not use this file except in compliance with the License.
|
||||
## You may obtain a copy of the License at
|
||||
##
|
||||
## http://www.apache.org/licenses/LICENSE-2.0
|
||||
##
|
||||
## Unless required by applicable law or agreed to in writing, software
|
||||
## distributed under the License is distributed on an "AS IS" BASIS,
|
||||
## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
||||
## See the License for the specific language governing permissions and
|
||||
## limitations under the License.
|
||||
|
||||
# \Phi(x)=2*(1-sqrt(1-x))/x
|
||||
@everywhere function Phi(x)
|
||||
if abs(x)>1e-5
|
||||
return 2*(1-sqrt(abs(1-x)))/x
|
||||
else
|
||||
return 1+x/4+x^2/8+5*x^3/64+7*x^4/128+21*x^5/512
|
||||
end
|
||||
end
|
||||
# \partial\Phi
|
||||
@everywhere function dPhi(x)
|
||||
#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))
|
||||
#end
|
||||
if abs(x)>1e-5
|
||||
return 1/(sqrt(abs(1-x))*x)*(1-x>=0 ? 1 : -1)-2*(1-sqrt(abs(1-x)))/x^2
|
||||
else
|
||||
return 1/4+x/4+15*x^2/64+7*x^3/32+105*x^4/512+99*x^5/512
|
||||
end
|
||||
end
|
||||
|
||||
# apply Phi to every element of a vector
|
||||
@everywhere function dotPhi(v)
|
||||
out=zeros(Float64,length(v))
|
||||
for i in 1:length(v)
|
||||
out[i]=Phi(v[i])
|
||||
end
|
||||
return out
|
||||
end
|
||||
@everywhere function dotdPhi(v)
|
||||
out=zeros(Float64,length(v))
|
||||
for i in 1:length(v)
|
||||
out[i]=dPhi(v[i])
|
||||
end
|
||||
return out
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user