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
|
/*
Copyright 2016 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.
*/
#include "zz_integral_double.h"
#include <math.h>
#define PI 3.1415926535897932385L
#define SQRT3 1.7320508075688772935L
// I for z1-z2
long double zz_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W){
long double xip,xiq,xiF,mp,mq,mF,zetap,zetaq,zetaF;
// alpha1 (store in xi)
xip=zz_alpha1_double(p1,p2);
xiq=zz_alpha1_double(q1,q2);
xiF=zz_alpha1_double(F1,F2);
// alpha2 (store in m)
mp=zz_alpha2_double(p1,p2);
mq=zz_alpha2_double(q1,q2);
mF=zz_alpha2_double(F1,F2);
zetap=2*t2*cosl(phi)*xip;
zetaq=2*t2*cosl(phi)*xiq;
zetaF=2*t2*cosl(phi)*xiF;
mp=W-2*t2*sinl(phi)*mp;
mq=W-2*t2*sinl(phi)*mq;
mF=W-2*t2*sinl(phi)*mF;
xip=sqrtl(mp*mp+2*t1*t1*xip);
xiq=sqrtl(mq*mq+2*t1*t1*xiq);
xiF=sqrtl(mF*mF+2*t1*t1*xiF);
return((xip+xiq+xiF)*(mp/xip+mq/xiq-mF/xiF-mp*mq*mF/xip/xiq/xiF)*(zetap+zetaq-zetaF)/(((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF))*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF)))/108);
}
// I for z1+z2
long double ZZ_I_double(long double p1, long double p2, long double q1, long double q2, long double F1, long double F2, long double t1, long double t2, long double phi, long double W){
long double xip,xiq,xiF,mp,mq,mF,zetap,zetaq,zetaF;
// alpha1 (store in xi)
xip=zz_alpha1_double(p1,p2);
xiq=zz_alpha1_double(q1,q2);
xiF=zz_alpha1_double(F1,F2);
// alpha2 (store in m)
mp=zz_alpha2_double(p1,p2);
mq=zz_alpha2_double(q1,q2);
mF=zz_alpha2_double(F1,F2);
zetap=2*t2*cosl(phi)*xip;
zetaq=2*t2*cosl(phi)*xiq;
zetaF=2*t2*cosl(phi)*xiF;
mp=W-2*t2*sinl(phi)*mp;
mq=W-2*t2*sinl(phi)*mq;
mF=W-2*t2*sinl(phi)*mF;
xip=sqrtl(mp*mp+2*t1*t1*xip);
xiq=sqrtl(mq*mq+2*t1*t1*xiq);
xiF=sqrtl(mF*mF+2*t1*t1*xiF);
return((1-mp/xip*mF/xiF-mq/xiq*mF/xiF+mp/xip*mq/xiq)*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)+(xip+xiq+xiF)*(xip+xiq+xiF))/(((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF))*((zetap+zetaq-zetaF)*(zetap+zetaq-zetaF)-(xip+xiq+xiF)*(xip+xiq+xiF)))/216);
}
// bar alpha_1(k1,k2)
long double zz_alpha1_double(long double k1, long double k2){
return(2*cosl(PI/3*k2)*(cosl(PI/SQRT3*k1)+cosl(PI/3*k2))+0.5);
}
// bar alpha_2(k1,k2)
long double zz_alpha2_double(long double k1, long double k2){
return(2*sinl(PI/3*k2)*(cosl(PI/SQRT3*k1)-cosl(PI/3*k2)));
}
|
