blob: 270519b955659b2d96c01fe914f131d8dc0fd575 (
plain)
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
|
o1 = dq**2
o2 = 1/o1
o3 = 1/dq
o4 = tan(qgrd(j))
o5 = 1/o4
o6 = deta**2
o7 = 1/o6
o8 = sin(qgrd(j))
o9 = o8**2
o10 = 1/o9
o11 = cos(qgrd(j))
o12 = o11**2
o13 = etagrd(i)**2
o14 = sigma**2
o15 = 1/o14
o16 = -(o13*o15)
o17 = -2.00000000000000d0*etagrd(i)*eta0*o15
o18 = eta0**2
o19 = -(o15*o18)
o20 = o16 + o17 + o19
o21 = exp(o20)
i22 = -2 + n
c The next lines causes a floating point exception on Alphas,
c if o22 is 0.0 but o8 is negative.
c This is because of y**x = exp(x*ln(y)).
c Changing the type of the exponent to integer fixes the problem.
c o22 = -2.00000000000000d0 + n
c o23 = o8**o22
o23 = o8**i22
o24 = n**2
o25 = 2.00000000000000d0*etagrd(i)*eta0*o15
o26 = o16 + o19 + o25
o27 = exp(o26)
o28 = o8**n
o29 = o14**2
o30 = 1/o29
o31 = o4**2
o32 = 1/o31
o33 = 5.0000000000000d-1*etagrd(i)
o34 = cosh(o33)
Cn(i,j) = o2 + 5.0000000000000d-1*o3*o5
Cs(i,j) = o2 - 5.0000000000000d-1*o3*o5
Cw(i,j) = o7
Cc(i,j) = -1.25000000000000d-1 - 1.25000000000000d-1*o10 - 2.000
& 00000000000d0*o2 - 2.50000000000000d-1*amp*n*o12*o21*o23 + 2.500
& 00000000000d-1*amp*o12*o21*o23*o24 - 2.50000000000000d-1*amp*n*o
& 12*o23*o27 + 2.50000000000000d-1*amp*o12*o23*o24*o27 - 2.5000000
& 0000000d-1*amp*n*o21*o28 - 5.0000000000000d-1*amp*o15*o21*o28 -
& 2.50000000000000d-1*amp*n*o27*o28 - 5.0000000000000d-1*amp*o15*o
& 27*o28 + 2.00000000000000d0*amp*etagrd(i)*eta0*o21*o28*o30 + amp*o13*o
& 21*o28*o30 + amp*o18*o21*o28*o30 - 2.00000000000000d0*amp*etagrd(i)*et
& a0*o27*o28*o30 + amp*o13*o27*o28*o30 + amp*o18*o27*o28*o30 + 1.2
& 5000000000000d-1*o32 - 2.00000000000000d0*o7
Ce(i,j) = o7
Rhs(i,j) = -2.50000000000000d-1*o34 + 2.50000000000000d-1*o10*o3
& 4 + 5.0000000000000d-1*amp*n*o12*o21*o23*o34 - 5.0000000000000d-
& 1*amp*o12*o21*o23*o24*o34 + 5.0000000000000d-1*amp*n*o12*o23*o27
& *o34 - 5.0000000000000d-1*amp*o12*o23*o24*o27*o34 + 5.0000000000
& 000d-1*amp*n*o21*o28*o34 + amp*o15*o21*o28*o34 + 5.0000000000000
& d-1*amp*n*o27*o28*o34 + amp*o15*o27*o28*o34 - 4.0000000000000d0*
& amp*etagrd(i)*eta0*o21*o28*o30*o34 - 2.00000000000000d0*amp*o13*o21*o2
& 8*o30*o34 - 2.00000000000000d0*amp*o18*o21*o28*o30*o34 + 4.00000
& 00000000d0*amp*etagrd(i)*eta0*o27*o28*o30*o34 - 2.00000000000000d0*amp
& *o13*o27*o28*o30*o34 - 2.00000000000000d0*amp*o18*o27*o28*o30*o3
& 4 - 2.50000000000000d-1*o32*o34
|
