This repository was archived by the owner on Jun 27, 2020. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathspecialfunctions.cpp
More file actions
141 lines (106 loc) · 11.1 KB
/
specialfunctions.cpp
File metadata and controls
141 lines (106 loc) · 11.1 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
/*
specialfunctions.cpp
Copyright (c) Michael Strickland
GNU General Public License (GPLv3)
See detailed text in license directory
*/
#include <cmath>
#include <iostream>
#include <fstream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <complex>
using namespace std;
#include "mpisolve.h"
#include "specialfunctions.h"
#define Power(a,b) pow((double)a,(double)b)
#define Log(a) log((double)a)
#define E exp(1)
// The Phi Function
inline double phiSmall(double r) {
return 6.356364414706291941e-25*Power(r,18)*(6.7858296469547001746e8 - 2.3279256e8*Log(r)) + 1.5134200987395933193e-26*Power(r,20)*(7.8066804966163335273e7 - 2.586584e7*Log(r)) + 7.3428721718687084502e-21*Power(r,16)*(1.7152136570933421512e7 - 6.12612e6*Log(r)) + 5.3476669453285429901e-15*Power(r,12)*(907957.56299608361835 - 360360.*Log(r)) + 2.970926080738079439e-17*Power(r,14)*(962011.56299608361835 - 360360.*Log(r)) + 9.0375571376052376533e-12*Power(r,10)*(64938.581768929509104 - 27720.*Log(r)) + 1.7496710618403740097e-8*Power(r,8)*(2679.7082622240685956 - 1260.*Log(r)) + 2.8344671201814058957e-6*Power(r,6)*(776.56942074135619855 - 420.*Log(r)) + 0.0011111111111111111111*Power(r,4)*(43.683530052954014182 - 30.*Log(r)) + 0.11111111111111111111*Power(r,2)*(2.2683530052954014182 - 3.*Log(r));
}
inline double phiMed(double r) {
return 0.96602365558836788831 + 0.0092523618649041681976*(-8. + r) - 0.0018940525872478075903*Power(-8. + r,2) + 0.00033240639398304988773*Power(-8. + r,3) - 0.000050871007785674343142*Power(-8. + r,4) + 6.7649512521654014838e-6*Power(-8. + r,5) - 7.7467140866580853702e-7*Power(-8. + r,6) + 7.4962910406047977489e-8*Power(-8. + r,7) - 5.8300876829970054706e-9*Power(-8. + r,8) + 3.0231708490081016772e-10*Power(-8. + r,9) + 3.1536456936207873081e-12*Power(-8. + r,10) - 3.5047790951879126093e-12*Power(-8. + r,11) + 6.0365878152222828005e-13*Power(-8. + r,12) - 7.6917232610971074521e-14*Power(-8. + r,13) + 8.5182961184006628674e-15*Power(-8. + r,14) - 8.7269091909431500681e-16*Power(-8. + r,15) + 8.5748841026634936066e-17*Power(-8. + r,16) - 8.2854177337578850044e-18*Power(-8. + r,17) + 8.013215406235426584e-19*Power(-8. + r,18) - 7.845320140747882741e-20*Power(-8. + r,19) + 7.8191161857969852618e-21*Power(-8. + r,20);
}
inline double phiLarge(double r) {
return 1. - 1.2804747411456e17/Power(r,20) - 3.76610217984e14/Power(r,18) - 1.3948526592e12/Power(r,16) - 6.7060224e9/Power(r,14) - 4.35456e7/Power(r,12) - 403200./Power(r,10) - 5760./Power(r,8) - 144./Power(r,6) - 8./Power(r,4) - 2./Power(r,2);
}
double Phi(double r) {
if (r <= 6) return phiSmall(r);
if (r > 6 && r < 13) return phiMed(r);
if (r>=13) return phiLarge(r);
return -1;
}
// The Psi1 Function
inline double Psi1_1_Small(double r){
return 0.5*r + 0.055555555555555555556*Power(r,3)*(-2.2683530052954014182 + 3.*Log(r)) + 0.00055555555555555555556*Power(r,5)*(-43.683530052954014182 + 30.*Log(r)) + 1.4172335600907029478e-6*Power(r,7)*(-776.56942074135619855 + 420.*Log(r)) + 8.7483553092018700484e-9*Power(r,9)*(-2679.7082622240685956 + 1260.*Log(r)) + 4.5187785688026188267e-12*Power(r,11)*(-64938.581768929509104 + 27720.*Log(r)) + 1.4854630403690397195e-17*Power(r,15)*(-962011.56299608361835 + 360360.*Log(r)) + 2.6738334726642714951e-15*Power(r,13)*(-907957.56299608361835 + 360360.*Log(r)) + 3.6714360859343542251e-21*Power(r,17)*(-1.7152136570933421512e7 + 6.12612e6*Log(r)) + 3.1781822073531459705e-25*Power(r,19)*(-6.7858296469547001746e8 + 2.3279256e8*Log(r));
}
inline double Psi1_1_Med(double r){
return 0.13590537764652844676 - 0.020021275253800616945*(-8. + r) + 0.0029500294165391462626*Power(-8. + r,2) - 0.00038259928230829575576*Power(-8. + r,3) + 0.0000372808341511724287*Power(-8. + r,4) - 1.6243011158244343644e-6*Power(-8. + r,5) - 2.8378999141946659382e-7*Power(-8. + r,6) + 8.7484062708712358552e-8*Power(-8. + r,7) - 1.4161104471035966862e-8*Power(-8. + r,8) + 1.7057755018952620644e-9*Power(-8. + r,9) - 1.6377312522488823309e-10*Power(-8. + r,10) + 1.2442293533941256783e-11*Power(-8. + r,11) - 6.6224557849495681557e-13*Power(-8. + r,12) + 5.8395396827701580568e-15*Power(-8. + r,13) + 4.3854318318828857907e-15*Power(-8. + r,14) - 7.6838438282307140647e-16*Power(-8. + r,15) + 9.3350095440617759139e-17*Power(-8. + r,16) - 9.7327495782859280157e-18*Power(-8. + r,17) + 9.3742270438477186861e-19*Power(-8. + r,18) - 8.6847964681856019558e-20*Power(-8. + r,19) + 7.9501359605514726576e-21*Power(-8. + r,20);
}
inline double Psi1_1_Large(double r){
return 6.402373705728e16/Power(r,19) + 1.88305108992e14/Power(r,17) + 6.974263296e11/Power(r,15) + 3.3530112e9/Power(r,13) + 2.17728e7/Power(r,11) + 201600./Power(r,9) + 2880./Power(r,7) + 72./Power(r,5) + 4./Power(r,3) + 1/r;
}
double psi1_1(double r) {
if (r <= 4) return Psi1_1_Small(r);
if (r > 4 && r < 14) return Psi1_1_Med(r);
if (r>=14) return Psi1_1_Large(r);
return -1;
}
inline double Psi1_2_Small(double r){
return -0.16666666666666666667 + 0.031870588947726682424*Power(r,2) + 0.0038070828840212936425*Power(r,4) + 0.00012963531232804592864*Power(r,6) + 2.2222840307597923849e-6*Power(r,8) + 2.3313730606766776239e-8*Power(r,10) + 1.6613089993632625133e-10*Power(r,12) + 8.591299661688259939e-13*Power(r,14) + 3.3766707701533263487e-15*Power(r,16) + 1.0437579040581881394e-17*Power(r,18) + 2.6054318816359345639e-20*Power(r,20) - 0.033333333333333333333*Power(r,2)*Log(r) - 0.0023809523809523809524*Power(r,4)*Log(r) - 0.00006613756613756613757*Power(r,6)*Log(r) - 1.002084335417668751e-6*Power(r,8)*Log(r) - 9.63542630209296876e-9*Power(r,10)*Log(r) - 6.42361753472864584e-11*Power(r,12)*Log(r) - 3.1488321248669832548e-13*Power(r,14)*Log(r) - 1.1837714755139034792e-15*Power(r,16)*Log(r) - 3.5231293914104270216e-18*Power(r,18)*Log(r) - 8.509974375387504883e-21*Power(r,20)*Log(r);
}
inline double Psi1_2_Med(double r){
return -0.0067413827273946599457 + 0.0010875998037452897687*(-12. + r) - 0.00012954793194127064687*Power(-12. + r,2) + 0.000013435401577769709693*Power(-12. + r,3) - 1.2707965149588749649e-6*Power(-12. + r,4) + 1.1129191549695214168e-7*Power(-12. + r,5) - 9.0455996105607973793e-9*Power(-12. + r,6) + 6.790200300786204485e-10*Power(-12. + r,7) - 4.6533092304693369172e-11*Power(-12. + r,8) + 2.8453161121616841998e-12*Power(-12. + r,9) - 1.4737612735202959067e-13*Power(-12. + r,10) + 5.47016964271159457e-15*Power(-12. + r,11) - 1.3638125698970240872e-18*Power(-12. + r,12) - 2.6199009725634488908e-17*Power(-12. + r,13) + 3.4409077531914510135e-18*Power(-12. + r,14) - 3.277190350118457649e-19*Power(-12. + r,15) + 2.6999082114796160642e-20*Power(-12. + r,16) - 2.0426080221425948919e-21*Power(-12. + r,17) + 1.4636193120195758959e-22*Power(-12. + r,18) - 1.0136545439995843953e-23*Power(-12. + r,19) + 6.8892595915062167481e-25*Power(-12. + r,20);
}
inline double Psi1_2_Large(double r){
return 3/(4.*Power(E,2*r)*Power(r,5)) + 3/Power(r,4) - 1/(2.*Power(E,2*r)*Power(r,3)) - Power(r,-2) - 1/(4.*Power(E,2*r)*Power(r,2));
}
double psi1_2(double r) {
if (r <= 3) return Psi1_2_Small(r);
if (r > 3 && r < 19) return Psi1_2_Med(r);
if (r>=19) return Psi1_2_Large(r);
return -1;
}
inline double Psi1(double r,double t) {
return 0.5 - 3*sin(t)*sin(t)*psi1_1(r)/2/r + 3*(1+3*cos(2*t))*psi1_2(r)/4;
}
// The Psi2 Function
inline double Psi2_1_Small(double r){
return 0.25*r - 0.041666666666666666667*Power(r,3) + 1.0569388174231281065e-18*Power(r,19)*(798.69116227224617869 - 280.*Log(r)) + 1.3382536530684678833e-13*Power(r,15)*(310.35017082786667278 - 120.*Log(r)) + 1.9680200780418645342e-15*Power(r,17)*(109.13633145242614583 - 40.*Log(r)) + 1.0036902398013509125e-10*Power(r,13)*(58.070034165573334556 - 24.*Log(r)) + 0.00055555555555555555556*Power(r,5)*(14.341765026477007091 - 15.*Log(r)) + 2.0876756987868098979e-8*Power(r,11)*(26.611940159709744201 - 12.*Log(r)) + 4.1335978835978835979e-6*Power(r,9)*(7.8403436896002177639 - 4.*Log(r)) + 0.00029761904761904761905*Power(r,7)*(3.1979496225778866597 - 2.*Log(r));
}
inline double Psi2_1_Med(double r){
return 0.008413723615770478 + 0.224967778471738063*r + 0.0369434258483233273*Power(r,2) - 0.0789499100348734325*Power(r,3) + 0.03271543254694100082*Power(r,4) - 0.0072819426449726113*Power(r,5) + 0.000887940602003804717*Power(r,6) - 5.4208645178527461e-6*Power(r,7) - 0.00002343699811431597517*Power(r,8) + 5.96255713043201882e-6*Power(r,9) - 9.40367493855617789e-7*Power(r,10) + 1.103982652930354787e-7*Power(r,11) - 1.0182758794827027921e-8*Power(r,12) + 7.517315207314647111e-10*Power(r,13) - 4.45167195175797198e-11*Power(r,14) + 2.095835211424293222e-12*Power(r,15) - 7.689144658891027814e-14*Power(r,16) + 2.1231926362574494762e-15*Power(r,17) - 4.156297438329627614e-17*Power(r,18) + 5.1466439416521259622e-19*Power(r,19) - 3.0329873366240464018e-21*Power(r,20);
}
inline double Psi2_1_Large(double r){
return 3.5213055381504e17/Power(r,19) + 9.4152554496e14/Power(r,17) + 3.1384184832e12/Power(r,15) + 1.34120448e10/Power(r,13) + 7.62048e7/Power(r,11) + 604800./Power(r,9) + 7200./Power(r,7) + 144./Power(r,5) + 6./Power(r,3) + 1/r;
}
double psi2_1(double r) {
if (r <= 5) return Psi2_1_Small(r);
if (r > 5 && r < 14) return Psi2_1_Med(r);
if (r>=14) return Psi2_1_Large(r);
return -1;
}
inline double Psi2_2_Small(double r){
return -0.083333333333333333333 + 0.0083333333333333333333*Power(r,2) - 0.0013083033467725515832*Power(r,4) - 0.00011310092079365439424*Power(r,6) - 3.0829049622852713896e-6*Power(r,8) - 4.4218604638010310288e-8*Power(r,10) - 3.9926820600399401372e-10*Power(r,12) - 2.4986690953848034005e-12*Power(r,14) - 1.1522404826658166351e-14*Power(r,16) - 4.0869533814474918822e-17*Power(r,18) - 1.1511694107977017916e-19*Power(r,20) + 3.8294884689243771973e-20*Power(r,4)* (3.10870812155904e16 + 1.7270600675328e15*Power(r,2) + 3.92513651712e13*Power(r,4) + 5.032226304e11*Power(r,6) + 4.19352192e9*Power(r,8) + 2.4667776e7*Power(r,10) + 108192.*Power(r,12) + 368.*Power(r,14) + Power(r,16))*Log(r);
}
inline double Psi2_2_Med(double r){
return -0.080454480144784501 - 0.00602738413377866263*r + 0.01481807672810240982*Power(r,2) - 0.005384527812250004745*Power(r,3) + 0.001116367841070016944*Power(r,4) - 0.0001524009859913323587*Power(r,5) + 0.00001309821798497723991*Power(r,6) - 3.48901092055444453e-7*Power(r,7) - 9.55173722279674947e-8*Power(r,8) + 1.947136662103399125e-8*Power(r,9) - 2.216686050000456054e-9*Power(r,10) + 1.836355223854977546e-10*Power(r,11) - 1.183326813340546626e-11*Power(r,12) + 6.063833665095452951e-13*Power(r,13) - 2.480127463641133029e-14*Power(r,14) + 8.030505474007715207e-16*Power(r,15) - 2.0189134732175212507e-17*Power(r,16) + 3.808098534812237147e-19*Power(r,17) - 5.078179382165779878e-21*Power(r,18) + 4.2733438967982267417e-23*Power(r,19) - 1.7078615998765512961e-25*Power(r,20);
}
inline double Psi2_2_Large(double r){
return 15/(16.*Power(E,2*r)*Power(r,5)) + 15/(4.*Power(r,4)) - 9/(16.*Power(E,2*r)*Power(r,3)) - Power(r,-2) - 5/(16.*Power(E,2*r)*Power(r,2)) - 1/(16.*Power(E,2*r)*r);
}
double psi2_2(double r) {
if (r <= 3) return Psi2_2_Small(r);
if (r > 3 && r < 19) return Psi2_2_Med(r);
if (r>=19) return Psi2_2_Large(r);
return -1;
}
double Psi2(double r,double t) {
return -1./3. + (2-6*cos(2*t))*psi2_1(r)/3/r - 2*(1+3*cos(2*t))*psi2_2(r);
}
double ImV(double r, double t, double xi) {
return -Phi(r) + xi*(Psi1(r,t)+Psi2(r,t));
}