Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module helm_polynomials
21 : use const_def, only: dp
22 :
23 : implicit none
24 :
25 : contains
26 :
27 : !..quintic hermite polynomial statement functions
28 : !..psi0 and its derivatives
29 182080 : pure real(dp) function psi0(z)
30 : real(dp), intent(in) :: z
31 182080 : psi0 = z*z*z* ( z * (-6.0d0*z + 15.0d0) - 10.0d0) + 1.0d0
32 182080 : end function psi0
33 :
34 182080 : pure real(dp) function dpsi0(z)
35 : real(dp), intent(in) :: z
36 182080 : dpsi0 = z*z* ( z * (-30.0d0*z + 60.0d0) - 30.0d0)
37 182080 : end function dpsi0
38 :
39 182080 : pure real(dp) function ddpsi0(z)
40 : real(dp), intent(in) :: z
41 182080 : ddpsi0 = z* ( z*( -120.0d0*z + 180.0d0) - 60.0d0)
42 182080 : end function ddpsi0
43 :
44 182080 : pure real(dp) function dddpsi0(z)
45 : real(dp), intent(in) :: z
46 182080 : dddpsi0 = z* ( -360.0d0*z + 360.0d0) - 60.0d0
47 182080 : end function dddpsi0
48 :
49 :
50 : !..psi1 and its derivatives
51 :
52 182080 : pure real(dp) function psi1(z)
53 : real(dp), intent(in) :: z
54 182080 : psi1 = z* (z*z * ( z * (-3.0d0*z + 8.0d0) - 6.0d0) + 1.0d0)
55 182080 : end function psi1
56 :
57 182080 : pure real(dp) function dpsi1(z)
58 : real(dp), intent(in) :: z
59 182080 : dpsi1 = z*z * ( z * (-15.0d0*z + 32.0d0) - 18.0d0) +1.0d0
60 182080 : end function dpsi1
61 :
62 182080 : pure real(dp) function ddpsi1(z)
63 : real(dp), intent(in) :: z
64 182080 : ddpsi1 = z * (z * (-60.0d0*z + 96.0d0) -36.0d0)
65 182080 : end function ddpsi1
66 :
67 182080 : pure real(dp) function dddpsi1(z)
68 : real(dp), intent(in) :: z
69 182080 : dddpsi1 = z * (-180.0d0*z + 192.0d0) - 36.0d0
70 182080 : end function dddpsi1
71 :
72 : !..psi2 and its derivatives
73 :
74 182080 : pure real(dp) function psi2(z)
75 : real(dp), intent(in) :: z
76 182080 : psi2 = 0.5d0*z*z*( z* ( z * (-z + 3.0d0) - 3.0d0) + 1.0d0)
77 182080 : end function psi2
78 :
79 182080 : pure real(dp) function dpsi2(z)
80 : real(dp), intent(in) :: z
81 182080 : dpsi2 = 0.5d0*z*( z*(z*(-5.0d0*z + 12.0d0) - 9.0d0) + 2.0d0)
82 182080 : end function dpsi2
83 :
84 182080 : pure real(dp) function ddpsi2(z)
85 : real(dp), intent(in) :: z
86 182080 : ddpsi2 = 0.5d0*(z*( z * (-20.0d0*z + 36.0d0) - 18.0d0) +2.0d0)
87 182080 : end function ddpsi2
88 :
89 182080 : pure real(dp) function dddpsi2(z)
90 : real(dp), intent(in) :: z
91 182080 : dddpsi2 = 0.5d0*(z * (-60.0d0*z + 72.0d0) - 18.0d0)
92 182080 : end function dddpsi2
93 :
94 : !..biquintic hermite polynomial statement function
95 :
96 455200 : pure real(dp) function h5(i, j, fi, &
97 : w0t, w1t, w2t, w0mt, w1mt, w2mt, &
98 : w0d, w1d, w2d, w0md, w1md, w2md)
99 : integer, intent(in) :: i, j
100 : real(dp), intent(in) :: fi(36)
101 : real(dp), intent(in) :: w0t, w1t, w2t, w0mt, w1mt, w2mt
102 : real(dp), intent(in) :: w0d, w1d, w2d, w0md, w1md, w2md
103 :
104 : h5 = fi(1) *w0d*w0t + fi(2) *w0md*w0t &
105 : + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt &
106 : + fi(5) *w0d*w1t + fi(6) *w0md*w1t &
107 : + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt &
108 : + fi(9) *w0d*w2t + fi(10) *w0md*w2t &
109 : + fi(11) *w0d*w2mt + fi(12) *w0md*w2mt &
110 : + fi(13) *w1d*w0t + fi(14) *w1md*w0t &
111 : + fi(15) *w1d*w0mt + fi(16) *w1md*w0mt &
112 : + fi(17) *w2d*w0t + fi(18) *w2md*w0t &
113 : + fi(19) *w2d*w0mt + fi(20) *w2md*w0mt &
114 : + fi(21) *w1d*w1t + fi(22) *w1md*w1t &
115 : + fi(23) *w1d*w1mt + fi(24) *w1md*w1mt &
116 : + fi(25) *w2d*w1t + fi(26) *w2md*w1t &
117 : + fi(27) *w2d*w1mt + fi(28) *w2md*w1mt &
118 : + fi(29) *w1d*w2t + fi(30) *w1md*w2t &
119 : + fi(31) *w1d*w2mt + fi(32) *w1md*w2mt &
120 : + fi(33) *w2d*w2t + fi(34) *w2md*w2t &
121 455200 : + fi(35) *w2d*w2mt + fi(36) *w2md*w2mt
122 455200 : end function h5
123 :
124 : !..cubic hermite polynomial statement functions
125 : !..psi0 & derivatives
126 :
127 182080 : pure real(dp) function xpsi0(z)
128 : real(dp), intent(in) :: z
129 182080 : xpsi0 = z * z * (2.0d0*z - 3.0d0) + 1.0d0
130 182080 : end function xpsi0
131 :
132 182080 : pure real(dp) function xdpsi0(z)
133 : real(dp), intent(in) :: z
134 182080 : xdpsi0 = z * (6.0d0*z - 6.0d0)
135 182080 : end function xdpsi0
136 :
137 182080 : pure real(dp) function xddpsi0(z)
138 : real(dp), intent(in) :: z
139 182080 : xddpsi0 = 12.0d0*z - 6.0d0
140 182080 : end function xddpsi0
141 :
142 182080 : pure real(dp) function xdddpsi0(z)
143 : real(dp), intent(in) :: z
144 182080 : xdddpsi0 = 12.0d0
145 182080 : end function xdddpsi0
146 :
147 : !..psi1 & derivatives
148 :
149 182080 : pure real(dp) function xpsi1(z)
150 : real(dp), intent(in) :: z
151 182080 : xpsi1 = z * ( z * (z - 2.0d0) + 1.0d0)
152 182080 : end function xpsi1
153 :
154 182080 : pure real(dp) function xdpsi1(z)
155 : real(dp), intent(in) :: z
156 182080 : xdpsi1 = z * (3.0d0*z - 4.0d0) + 1.0d0
157 182080 : end function xdpsi1
158 :
159 182080 : pure real(dp) function xddpsi1(z)
160 : real(dp), intent(in) :: z
161 182080 : xddpsi1 = 6.0d0*z - 4.0d0
162 182080 : end function xddpsi1
163 :
164 182080 : pure real(dp) function xdddpsi1(z)
165 : real(dp), intent(in) :: z
166 182080 : xdddpsi1 = 6.0d0
167 182080 : end function xdddpsi1
168 :
169 : !..bicubic hermite polynomial statement function
170 910400 : pure real(dp) function h3(i, j, fi, &
171 : w0t, w1t, w0mt, w1mt, w0d, w1d, w0md, w1md)
172 : integer, intent(in) :: i, j
173 : real(dp), intent(in) :: fi(36)
174 : real(dp), intent(in) :: w0t, w1t, w0mt, w1mt, w0d, w1d, w0md, w1md
175 : h3 = fi(1) *w0d*w0t + fi(2) *w0md*w0t &
176 : + fi(3) *w0d*w0mt + fi(4) *w0md*w0mt &
177 : + fi(5) *w0d*w1t + fi(6) *w0md*w1t &
178 : + fi(7) *w0d*w1mt + fi(8) *w0md*w1mt &
179 : + fi(9) *w1d*w0t + fi(10) *w1md*w0t &
180 : + fi(11) *w1d*w0mt + fi(12) *w1md*w0mt &
181 : + fi(13) *w1d*w1t + fi(14) *w1md*w1t &
182 910400 : + fi(15) *w1d*w1mt + fi(16) *w1md*w1mt
183 910400 : end function h3
184 :
185 : end module helm_polynomials
|
