1088 integer,
intent(in) :: weno_dir
1089 type(int_bounds_info),
intent(in) :: is
1091 real(wp),
pointer,
dimension(:) :: s_cb => null()
1092 type(int_bounds_info) :: bc_s
1100 if (weno_dir == 1)
then
1101 s = m; s_cb => x_cb; bc_s = bc_x
1102 else if (weno_dir == 2)
then
1103 s = n; s_cb => y_cb; bc_s = bc_y
1105 s = p; s_cb => z_cb; bc_s = bc_z
1108# 192 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
1110 if (weno_dir == 1)
then
1111 if (weno_order == 3)
then
1112 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
1114 poly_coef_cbr_x(i + 1, 0, 0) = (s_cb(i) - s_cb(i + 1))/(s_cb(i) - s_cb(i + 2))
1115 poly_coef_cbr_x(i + 1, 1, 0) = (s_cb(i) - s_cb(i + 1))/(s_cb(i - 1) - s_cb(i + 1))
1121 d_cbr_x(0, i + 1) = (s_cb(i - 1) - s_cb(i + 1))/(s_cb(i - 1) - s_cb(i + 2))
1122 d_cbl_x(0, i + 1) = (s_cb(i - 1) - s_cb(i))/(s_cb(i - 1) - s_cb(i + 2))
1128 beta_coef_x(i + 1, 0, 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp/(s_cb(i) - s_cb(i + 2))**2._wp
1129 beta_coef_x(i + 1, 1, 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp/(s_cb(i - 1) - s_cb(i + 1))**2._wp
1134 if (null_weights)
then
1135 if (bc_s%beg == bc_riemann_extrap)
then
1140 if (bc_s%end == bc_riemann_extrap)
then
1148 else if (weno_order == 5)
then
1149 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
1152 & 0) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i + 2)))/((s_cb(i) - s_cb(i &
1153 & + 3))*(s_cb(i + 3) - s_cb(i + 1)))
1155 & 0) = ((s_cb(i - 1) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i)))/((s_cb(i - 1) &
1156 & - s_cb(i + 2))*(s_cb(i + 2) - s_cb(i)))
1158 & 1) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i + 2)))/((s_cb(i - 1) &
1159 & - s_cb(i + 1))*(s_cb(i - 1) - s_cb(i + 2)))
1161 & 1) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))/((s_cb(i - 2) &
1162 & - s_cb(i))*(s_cb(i - 2) - s_cb(i + 1)))
1164 & 0) = ((s_cb(i + 1) - s_cb(i))*(s_cb(i) - s_cb(i + 2)))/((s_cb(i) - s_cb(i + 3)) &
1165 & *(s_cb(i + 3) - s_cb(i + 1)))
1167 & 0) = ((s_cb(i) - s_cb(i - 1))*(s_cb(i) - s_cb(i + 1)))/((s_cb(i - 1) - s_cb(i &
1168 & + 2))*(s_cb(i) - s_cb(i + 2)))
1170 & 1) = ((s_cb(i + 1) - s_cb(i))*(s_cb(i) - s_cb(i + 2)))/((s_cb(i - 1) - s_cb(i &
1171 & + 1))*(s_cb(i - 1) - s_cb(i + 2)))
1173 & 1) = ((s_cb(i - 1) - s_cb(i))*(s_cb(i) - s_cb(i + 1)))/((s_cb(i - 2) - s_cb(i)) &
1174 & *(s_cb(i - 2) - s_cb(i + 1)))
1177 & 1) = ((s_cb(i) - s_cb(i + 2)) + (s_cb(i + 1) - s_cb(i + 3)))/((s_cb(i) - s_cb(i &
1178 & + 2))*(s_cb(i) - s_cb(i + 3)))*((s_cb(i) - s_cb(i + 1)))
1180 & 0) = ((s_cb(i - 2) - s_cb(i + 1)) + (s_cb(i - 1) - s_cb(i + 1)))/((s_cb(i - 1) &
1181 & - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 2)))*((s_cb(i + 1) - s_cb(i)))
1183 & 1) = ((s_cb(i) - s_cb(i + 2)) + (s_cb(i) - s_cb(i + 3)))/((s_cb(i) - s_cb(i + 2)) &
1184 & *(s_cb(i) - s_cb(i + 3)))*((s_cb(i + 1) - s_cb(i)))
1186 & 0) = ((s_cb(i - 2) - s_cb(i)) + (s_cb(i - 1) - s_cb(i + 1)))/((s_cb(i - 2) &
1187 & - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))*((s_cb(i) - s_cb(i + 1)))
1191 & i + 1) = ((s_cb(i - 2) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))/((s_cb(i - 2) &
1192 & - s_cb(i + 3))*(s_cb(i + 3) - s_cb(i - 1)))
1194 & i + 1) = ((s_cb(i + 1) - s_cb(i + 2))*(s_cb(i + 1) - s_cb(i + 3)))/((s_cb(i - 2) &
1195 & - s_cb(i + 2))*(s_cb(i - 2) - s_cb(i + 3)))
1197 & i + 1) = ((s_cb(i - 2) - s_cb(i))*(s_cb(i) - s_cb(i - 1)))/((s_cb(i - 2) - s_cb(i + 3)) &
1198 & *(s_cb(i + 3) - s_cb(i - 1)))
1200 & i + 1) = ((s_cb(i) - s_cb(i + 2))*(s_cb(i) - s_cb(i + 3)))/((s_cb(i - 2) - s_cb(i + 2)) &
1201 & *(s_cb(i - 2) - s_cb(i + 3)))
1208 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1209 & + (s_cb(i + 1) - s_cb(i))*(s_cb(i + 2) - s_cb(i + 1)) + (s_cb(i + 2) - s_cb(i + 1)) &
1210 & **2._wp)/((s_cb(i) - s_cb(i + 3))**2._wp*(s_cb(i + 1) - s_cb(i + 3))**2._wp)
1213 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(19._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1214 & - (s_cb(i + 1) - s_cb(i))*(s_cb(i + 3) - s_cb(i + 1)) + 2._wp*(s_cb(i + 2) - s_cb(i)) &
1215 & *((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1))))/((s_cb(i) - s_cb(i + 2)) &
1216 & *(s_cb(i) - s_cb(i + 3))**2._wp*(s_cb(i + 3) - s_cb(i + 1)))
1219 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1220 & + (s_cb(i + 1) - s_cb(i))*((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1))) &
1221 & + ((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1)))**2._wp)/((s_cb(i) - s_cb(i &
1222 & + 2))**2._wp*(s_cb(i) - s_cb(i + 3))**2._wp)
1225 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1226 & + (s_cb(i) - s_cb(i - 1))**2._wp + (s_cb(i) - s_cb(i - 1))*(s_cb(i + 1) - s_cb(i))) &
1227 & /((s_cb(i - 1) - s_cb(i + 2))**2._wp*(s_cb(i) - s_cb(i + 2))**2._wp)
1230 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*((s_cb(i) - s_cb(i + 1))*((s_cb(i) &
1231 & - s_cb(i - 1)) + 20._wp*(s_cb(i + 1) - s_cb(i))) + (2._wp*(s_cb(i) - s_cb(i - 1)) &
1232 & + (s_cb(i + 1) - s_cb(i)))*(s_cb(i + 2) - s_cb(i)))/((s_cb(i + 1) - s_cb(i - 1)) &
1233 & *(s_cb(i - 1) - s_cb(i + 2))**2._wp*(s_cb(i + 2) - s_cb(i)))
1236 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1237 & + (s_cb(i + 1) - s_cb(i))*(s_cb(i + 2) - s_cb(i + 1)) + (s_cb(i + 2) - s_cb(i + 1)) &
1238 & **2._wp)/((s_cb(i - 1) - s_cb(i + 1))**2._wp*(s_cb(i - 1) - s_cb(i + 2))**2._wp)
1241 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(12._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1242 & + ((s_cb(i) - s_cb(i - 2)) + (s_cb(i) - s_cb(i - 1)))**2._wp + 3._wp*((s_cb(i) &
1243 & - s_cb(i - 2)) + (s_cb(i) - s_cb(i - 1)))*(s_cb(i + 1) - s_cb(i)))/((s_cb(i - 2) &
1244 & - s_cb(i + 1))**2._wp*(s_cb(i - 1) - s_cb(i + 1))**2._wp)
1247 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(19._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1248 & + ((s_cb(i) - s_cb(i - 2))*(s_cb(i) - s_cb(i + 1))) + 2._wp*(s_cb(i + 1) - s_cb(i &
1249 & - 1))*((s_cb(i) - s_cb(i - 2)) + (s_cb(i + 1) - s_cb(i - 1))))/((s_cb(i - 2) &
1250 & - s_cb(i))*(s_cb(i - 2) - s_cb(i + 1))**2._wp*(s_cb(i + 1) - s_cb(i - 1)))
1253 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1254 & + (s_cb(i) - s_cb(i - 1))**2._wp + (s_cb(i) - s_cb(i - 1))*(s_cb(i + 1) - s_cb(i))) &
1255 & /((s_cb(i - 2) - s_cb(i))**2._wp*(s_cb(i - 2) - s_cb(i + 1))**2._wp)
1260 if (null_weights)
then
1261 if (bc_s%beg == bc_riemann_extrap)
then
1268 if (bc_s%end == bc_riemann_extrap)
then
1270 & s - 1)/sum(
d_cbr_x(:,s - 1))
1272 & s - 1)/sum(
d_cbl_x(:,s - 1))
1278 if (.not. teno)
then
1279 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
1292 w = s_cb(i - 3:i + 4) - s_cb(i)
1294 & i + 1) = ((w(5) - w(6))*(w(5) - w(7))*(w(5) - w(8)))/((w(1) - w(6))*(w(1) - w(7)) &
1297 & i + 1) = ((w(1) - w(5))*(w(5) - w(7))*(w(5) - w(8))*(w(1)*w(2) - w(1)*w(6) - w(1) &
1298 & *w(7) - w(2)*w(6) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) + w(6)*w(7) + w(6)*w(8) + w(7) &
1299 & *w(8) + w(1)**2 + w(2)**2))/((w(1) - w(6))*(w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7)) &
1302 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(5) - w(8))*(w(1)*w(2) + w(1)*w(3) + w(2) &
1303 & *w(3) - w(1)*w(7) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) - w(3)*w(7) - w(3)*w(8) + w(7) &
1304 & *w(8) + w(7)**2 + w(8)**2))/((w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7))*(w(2) - w(8)) &
1307 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(3) - w(5)))/((w(1) - w(8))*(w(2) - w(8)) &
1314 w(1) = s_cb(i + 4) - s_cb(i)
1315 w(2) = s_cb(i + 3) - s_cb(i)
1316 w(3) = s_cb(i + 2) - s_cb(i)
1317 w(4) = s_cb(i + 1) - s_cb(i)
1318 w(5) = s_cb(i) - s_cb(i)
1319 w(6) = s_cb(i - 1) - s_cb(i)
1320 w(7) = s_cb(i - 2) - s_cb(i)
1321 w(8) = s_cb(i - 3) - s_cb(i)
1323 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(3) - w(5)))/((w(1) - w(8))*(w(2) - w(8)) &
1326 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(5) - w(8))*(w(1)*w(2) + w(1)*w(3) + w(2) &
1327 & *w(3) - w(1)*w(7) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) - w(3)*w(7) - w(3)*w(8) + w(7) &
1328 & *w(8) + w(7)**2 + w(8)**2))/((w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7))*(w(2) - w(8)) &
1331 & i + 1) = ((w(1) - w(5))*(w(5) - w(7))*(w(5) - w(8))*(w(1)*w(2) - w(1)*w(6) - w(1) &
1332 & *w(7) - w(2)*w(6) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) + w(6)*w(7) + w(6)*w(8) + w(7) &
1333 & *w(8) + w(1)**2 + w(2)**2))/((w(1) - w(6))*(w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7)) &
1336 & i + 1) = ((w(5) - w(6))*(w(5) - w(7))*(w(5) - w(8)))/((w(1) - w(6))*(w(1) - w(7)) &
1340 y = s_cb(i + 1:i + 4) - s_cb(i:i + 3)
1342 & 0) = (y(1)*y(2)*(y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
1343 & + y(2) + y(3) + y(4)))
1345 & 1) = -(y(1)*y(2)*(3*y(2)**2 + 6*y(2)*y(3) + 3*y(2)*y(4) + 2*y(1)*y(2) &
1346 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
1347 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1349 & 2) = (y(1)*(y(1)**2 + 3*y(1)*y(2) + 2*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
1350 & + 4*y(2)*y(3) + 2*y(4)*y(2) + y(3)**2 + y(4)*y(3)))/((y(1) + y(2))*(y(1) &
1351 & + y(2) + y(3))*(y(1) + y(2) + y(3) + y(4)))
1353 y = s_cb(i:i + 3) - s_cb(i - 1:i + 2)
1355 & 0) = -(y(2)*y(3)*(y(1) + y(2)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
1356 & + y(2) + y(3) + y(4)))
1358 & 1) = (y(2)*(y(1) + y(2))*(y(2)**2 + 4*y(2)*y(3) + 2*y(2)*y(4) + y(1)*y(2) &
1359 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
1360 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1362 & 2) = (y(2)*y(3)*(y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
1363 & + y(2) + y(3) + y(4)))
1365 y = s_cb(i - 1:i + 2) - s_cb(i - 2:i + 1)
1367 & 0) = (y(3)*(y(2) + y(3))*(y(1) + y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) &
1368 & + y(4))*(y(1) + y(2) + y(3) + y(4)))
1370 & 1) = (y(3)*y(4)*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
1371 & + 6*y(2)*y(3) + 2*y(4)*y(2) + 3*y(3)**2 + 2*y(4)*y(3)))/((y(2) + y(3))*(y(1) &
1372 & + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1374 & 2) = -(y(3)*y(4)*(y(2) + y(3)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
1375 & + y(2) + y(3) + y(4)))
1377 y = s_cb(i - 2:i + 1) - s_cb(i - 3:i)
1379 & 0) = (y(4)*(y(2)**2 + 4*y(2)*y(3) + 4*y(2)*y(4) + y(1)*y(2) + 3*y(3)**2 &
1380 & + 6*y(3)*y(4) + 2*y(1)*y(3) + 3*y(4)**2 + 2*y(1)*y(4)))/((y(3) + y(4))*(y(2) &
1381 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1383 & 1) = -(y(4)*(y(3) + y(4))*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + 2*y(1)*y(4) &
1384 & + 3*y(2)**2 + 6*y(2)*y(3) + 4*y(2)*y(4) + 3*y(3)**2 + 4*y(3)*y(4) + y(4)**2)) &
1385 & /((y(2) + y(3))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
1388 & 2) = (y(4)*(y(3) + y(4))*(y(2) + y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) &
1389 & + y(3))*(y(1) + y(2) + y(3) + y(4)))
1392 y(1) = s_cb(i + 1) - s_cb(i)
1393 y(2) = s_cb(i) - s_cb(i - 1)
1394 y(3) = s_cb(i - 1) - s_cb(i - 2)
1395 y(4) = s_cb(i - 2) - s_cb(i - 3)
1397 & 2) = (y(1)*y(2)*(y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
1398 & + y(2) + y(3) + y(4)))
1400 & 1) = -(y(1)*y(2)*(3*y(2)**2 + 6*y(2)*y(3) + 3*y(2)*y(4) + 2*y(1)*y(2) &
1401 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
1402 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1404 & 0) = (y(1)*(y(1)**2 + 3*y(1)*y(2) + 2*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
1405 & + 4*y(2)*y(3) + 2*y(4)*y(2) + y(3)**2 + y(4)*y(3)))/((y(1) + y(2))*(y(1) &
1406 & + y(2) + y(3))*(y(1) + y(2) + y(3) + y(4)))
1409 y(1) = s_cb(i + 2) - s_cb(i + 1)
1410 y(2) = s_cb(i + 1) - s_cb(i)
1411 y(3) = s_cb(i) - s_cb(i - 1)
1412 y(4) = s_cb(i - 1) - s_cb(i - 2)
1414 & 2) = -(y(2)*y(3)*(y(1) + y(2)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
1415 & + y(2) + y(3) + y(4)))
1417 & 1) = (y(2)*(y(1) + y(2))*(y(2)**2 + 4*y(2)*y(3) + 2*y(2)*y(4) + y(1)*y(2) &
1418 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
1419 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1421 & 0) = (y(2)*y(3)*(y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
1422 & + y(2) + y(3) + y(4)))
1425 y(1) = s_cb(i + 3) - s_cb(i + 2)
1426 y(2) = s_cb(i + 2) - s_cb(i + 1)
1427 y(3) = s_cb(i + 1) - s_cb(i)
1428 y(4) = s_cb(i) - s_cb(i - 1)
1430 & 2) = (y(3)*(y(2) + y(3))*(y(1) + y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) &
1431 & + y(4))*(y(1) + y(2) + y(3) + y(4)))
1433 & 1) = (y(3)*y(4)*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
1434 & + 6*y(2)*y(3) + 2*y(4)*y(2) + 3*y(3)**2 + 2*y(4)*y(3)))/((y(2) + y(3))*(y(1) &
1435 & + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1437 & 0) = -(y(3)*y(4)*(y(2) + y(3)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
1438 & + y(2) + y(3) + y(4)))
1441 y(1) = s_cb(i + 4) - s_cb(i + 3)
1442 y(2) = s_cb(i + 3) - s_cb(i + 2)
1443 y(3) = s_cb(i + 2) - s_cb(i + 1)
1444 y(4) = s_cb(i + 1) - s_cb(i)
1446 & 2) = (y(4)*(y(2)**2 + 4*y(2)*y(3) + 4*y(2)*y(4) + y(1)*y(2) + 3*y(3)**2 &
1447 & + 6*y(3)*y(4) + 2*y(1)*y(3) + 3*y(4)**2 + 2*y(1)*y(4)))/((y(3) + y(4))*(y(2) &
1448 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
1450 & 1) = -(y(4)*(y(3) + y(4))*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + 2*y(1)*y(4) &
1451 & + 3*y(2)**2 + 6*y(2)*y(3) + 4*y(2)*y(4) + 3*y(3)**2 + 4*y(3)*y(4) + y(4)**2)) &
1452 & /((y(2) + y(3))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
1455 & 0) = (y(4)*(y(3) + y(4))*(y(2) + y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) &
1456 & + y(3))*(y(1) + y(2) + y(3) + y(4)))
1461 y = s_cb(i - 2:i + 1) - s_cb(i - 3:i)
1463 & 0) = (4*y(4)**2*(5*y(1)**2*y(2)**2 + 20*y(1)**2*y(2)*y(3) + 15*y(1)**2*y(2)*y(4) &
1464 & + 20*y(1)**2*y(3)**2 + 30*y(1)**2*y(3)*y(4) + 60*y(1)**2*y(4)**2 + 10*y(1)*y(2) &
1465 & **3 + 60*y(1)*y(2)**2*y(3) + 45*y(1)*y(2)**2*y(4) + 110*y(1)*y(2)*y(3)**2 &
1466 & + 165*y(1)*y(2)*y(3)*y(4) + 260*y(1)*y(2)*y(4)**2 + 60*y(1)*y(3)**3 + 135*y(1) &
1467 & *y(3)**2*y(4) + 400*y(1)*y(3)*y(4)**2 + 225*y(1)*y(4)**3 + 5*y(2)**4 + 40*y(2) &
1468 & **3*y(3) + 30*y(2)**3*y(4) + 110*y(2)**2*y(3)**2 + 165*y(2)**2*y(3)*y(4) &
1469 & + 260*y(2)**2*y(4)**2 + 120*y(2)*y(3)**3 + 270*y(2)*y(3)**2*y(4) + 800*y(2)*y(3) &
1470 & *y(4)**2 + 450*y(2)*y(4)**3 + 45*y(3)**4 + 135*y(3)**3*y(4) + 600*y(3)**2*y(4) &
1471 & **2 + 675*y(3)*y(4)**3 + 996*y(4)**4))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4)) &
1472 & **2*(y(1) + y(2) + y(3) + y(4))**2)
1474 & 1) = -(4*y(4)**2*(10*y(1)**3*y(2)*y(3) + 5*y(1)**3*y(2)*y(4) + 20*y(1)**3*y(3) &
1475 & **2 + 25*y(1)**3*y(3)*y(4) + 105*y(1)**3*y(4)**2 + 40*y(1)**2*y(2)**2*y(3) &
1476 & + 20*y(1)**2*y(2)**2*y(4) + 130*y(1)**2*y(2)*y(3)**2 + 155*y(1)**2*y(2)*y(3)*y(4) &
1477 & + 535*y(1)**2*y(2)*y(4)**2 + 90*y(1)**2*y(3)**3 + 165*y(1)**2*y(3)**2*y(4) &
1478 & + 790*y(1)**2*y(3)*y(4)**2 + 415*y(1)**2*y(4)**3 + 60*y(1)*y(2)**3*y(3) + 30*y(1) &
1479 & *y(2)**3*y(4) + 270*y(1)*y(2)**2*y(3)**2 + 315*y(1)*y(2)**2*y(3)*y(4) + 975*y(1) &
1480 & *y(2)**2*y(4)**2 + 360*y(1)*y(2)*y(3)**3 + 645*y(1)*y(2)*y(3)**2*y(4) + 2850*y(1) &
1481 & *y(2)*y(3)*y(4)**2 + 1460*y(1)*y(2)*y(4)**3 + 150*y(1)*y(3)**4 + 360*y(1)*y(3) &
1482 & **3*y(4) + 2000*y(1)*y(3)**2*y(4)**2 + 2005*y(1)*y(3)*y(4)**3 + 2077*y(1)*y(4) &
1483 & **4 + 30*y(2)**4*y(3) + 15*y(2)**4*y(4) + 180*y(2)**3*y(3)**2 + 210*y(2)**3*y(3) &
1484 & *y(4) + 650*y(2)**3*y(4)**2 + 360*y(2)**2*y(3)**3 + 645*y(2)**2*y(3)**2*y(4) &
1485 & + 2850*y(2)**2*y(3)*y(4)**2 + 1460*y(2)**2*y(4)**3 + 300*y(2)*y(3)**4 + 720*y(2) &
1486 & *y(3)**3*y(4) + 4000*y(2)*y(3)**2*y(4)**2 + 4010*y(2)*y(3)*y(4)**3 + 4154*y(2) &
1487 & *y(4)**4 + 90*y(3)**5 + 270*y(3)**4*y(4) + 1800*y(3)**3*y(4)**2 + 2655*y(3) &
1488 & **2*y(4)**3 + 4464*y(3)*y(4)**4 + 1767*y(4)**5))/(5*(y(2) + y(3))*(y(3) + y(4)) &
1489 & *(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
1491 & 2) = (4*y(4)**2*(10*y(2)**3*y(3) + 5*y(2)**3*y(4) + 50*y(2)**2*y(3)**2 + 60*y(2) &
1492 & **2*y(3)*y(4) + 10*y(1)*y(2)**2*y(3) + 215*y(2)**2*y(4)**2 + 5*y(1)*y(2)**2*y(4) &
1493 & + 70*y(2)*y(3)**3 + 130*y(2)*y(3)**2*y(4) + 30*y(1)*y(2)*y(3)**2 + 775*y(2)*y(3) &
1494 & *y(4)**2 + 35*y(1)*y(2)*y(3)*y(4) + 415*y(2)*y(4)**3 + 110*y(1)*y(2)*y(4)**2 &
1495 & + 30*y(3)**4 + 75*y(3)**3*y(4) + 20*y(1)*y(3)**3 + 665*y(3)**2*y(4)**2 + 35*y(1) &
1496 & *y(3)**2*y(4) + 725*y(3)*y(4)**3 + 220*y(1)*y(3)*y(4)**2 + 1767*y(4)**4 &
1497 & + 105*y(1)*y(4)**3))/(5*(y(1) + y(2))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) &
1498 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
1500 & 3) = (4*y(4)**2*(5*y(1)**4*y(3)**2 + 5*y(1)**4*y(3)*y(4) + 50*y(1)**4*y(4)**2 &
1501 & + 30*y(1)**3*y(2)*y(3)**2 + 30*y(1)**3*y(2)*y(3)*y(4) + 300*y(1)**3*y(2)*y(4)**2 &
1502 & + 30*y(1)**3*y(3)**3 + 45*y(1)**3*y(3)**2*y(4) + 415*y(1)**3*y(3)*y(4)**2 &
1503 & + 200*y(1)**3*y(4)**3 + 75*y(1)**2*y(2)**2*y(3)**2 + 75*y(1)**2*y(2)**2*y(3)*y(4) &
1504 & + 750*y(1)**2*y(2)**2*y(4)**2 + 150*y(1)**2*y(2)*y(3)**3 + 225*y(1)**2*y(2)*y(3) &
1505 & **2*y(4) + 2075*y(1)**2*y(2)*y(3)*y(4)**2 + 1000*y(1)**2*y(2)*y(4)**3 + 75*y(1) &
1506 & **2*y(3)**4 + 150*y(1)**2*y(3)**3*y(4) + 1390*y(1)**2*y(3)**2*y(4)**2 + 1315*y(1) &
1507 & **2*y(3)*y(4)**3 + 1081*y(1)**2*y(4)**4 + 90*y(1)*y(2)**3*y(3)**2 + 90*y(1)*y(2) &
1508 & **3*y(3)*y(4) + 900*y(1)*y(2)**3*y(4)**2 + 270*y(1)*y(2)**2*y(3)**3 + 405*y(1) &
1509 & *y(2)**2*y(3)**2*y(4) + 3735*y(1)*y(2)**2*y(3)*y(4)**2 + 1800*y(1)*y(2)**2*y(4) &
1510 & **3 + 270*y(1)*y(2)*y(3)**4 + 540*y(1)*y(2)*y(3)**3*y(4) + 5025*y(1)*y(2)*y(3) &
1511 & **2*y(4)**2 + 4755*y(1)*y(2)*y(3)*y(4)**3 + 4224*y(1)*y(2)*y(4)**4 + 90*y(1)*y(3) &
1512 & **5 + 225*y(1)*y(3)**4*y(4) + 2190*y(1)*y(3)**3*y(4)**2 + 3060*y(1)*y(3)**2*y(4) &
1513 & **3 + 4529*y(1)*y(3)*y(4)**4 + 1762*y(1)*y(4)**5 + 45*y(2)**4*y(3)**2 + 45*y(2) &
1514 & **4*y(3)*y(4) + 450*y(2)**4*y(4)**2 + 180*y(2)**3*y(3)**3 + 270*y(2)**3*y(3) &
1515 & **2*y(4) + 2490*y(2)**3*y(3)*y(4)**2 + 1200*y(2)**3*y(4)**3 + 270*y(2)**2*y(3) &
1516 & **4 + 540*y(2)**2*y(3)**3*y(4) + 5025*y(2)**2*y(3)**2*y(4)**2 + 4755*y(2)**2*y(3) &
1517 & *y(4)**3 + 4224*y(2)**2*y(4)**4 + 180*y(2)*y(3)**5 + 450*y(2)*y(3)**4*y(4) &
1518 & + 4380*y(2)*y(3)**3*y(4)**2 + 6120*y(2)*y(3)**2*y(4)**3 + 9058*y(2)*y(3)*y(4)**4 &
1519 & + 3524*y(2)*y(4)**5 + 45*y(3)**6 + 135*y(3)**5*y(4) + 1395*y(3)**4*y(4)**2 &
1520 & + 2565*y(3)**3*y(4)**3 + 4884*y(3)**2*y(4)**4 + 3624*y(3)*y(4)**5 + 831*y(4)**6)) &
1521 & /(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
1522 & + y(3) + y(4))**2)
1524 & 4) = -(4*y(4)**2*(10*y(1)**2*y(2)*y(3)**2 + 10*y(1)**2*y(2)*y(3)*y(4) + 100*y(1) &
1525 & **2*y(2)*y(4)**2 + 10*y(1)**2*y(3)**3 + 15*y(1)**2*y(3)**2*y(4) + 205*y(1) &
1526 & **2*y(3)*y(4)**2 + 100*y(1)**2*y(4)**3 + 30*y(1)*y(2)**2*y(3)**2 + 30*y(1)*y(2) &
1527 & **2*y(3)*y(4) + 300*y(1)*y(2)**2*y(4)**2 + 60*y(1)*y(2)*y(3)**3 + 90*y(1)*y(2) &
1528 & *y(3)**2*y(4) + 1030*y(1)*y(2)*y(3)*y(4)**2 + 500*y(1)*y(2)*y(4)**3 + 30*y(1) &
1529 & *y(3)**4 + 60*y(1)*y(3)**3*y(4) + 835*y(1)*y(3)**2*y(4)**2 + 805*y(1)*y(3)*y(4) &
1530 & **3 + 1762*y(1)*y(4)**4 + 30*y(2)**3*y(3)**2 + 30*y(2)**3*y(3)*y(4) + 300*y(2) &
1531 & **3*y(4)**2 + 90*y(2)**2*y(3)**3 + 135*y(2)**2*y(3)**2*y(4) + 1445*y(2)**2*y(3) &
1532 & *y(4)**2 + 700*y(2)**2*y(4)**3 + 90*y(2)*y(3)**4 + 180*y(2)*y(3)**3*y(4) &
1533 & + 2205*y(2)*y(3)**2*y(4)**2 + 2115*y(2)*y(3)*y(4)**3 + 3624*y(2)*y(4)**4 &
1534 & + 30*y(3)**5 + 75*y(3)**4*y(4) + 1060*y(3)**3*y(4)**2 + 1515*y(3)**2*y(4)**3 &
1535 & + 3824*y(3)*y(4)**4 + 1662*y(4)**5))/(5*(y(1) + y(2))*(y(2) + y(3))*(y(1) + y(2) &
1536 & + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
1538 & 5) = (4*y(4)**2*(5*y(2)**2*y(3)**2 + 5*y(2)**2*y(3)*y(4) + 50*y(2)**2*y(4)**2 &
1539 & + 10*y(2)*y(3)**3 + 15*y(2)*y(3)**2*y(4) + 205*y(2)*y(3)*y(4)**2 + 100*y(2)*y(4) &
1540 & **3 + 5*y(3)**4 + 10*y(3)**3*y(4) + 205*y(3)**2*y(4)**2 + 200*y(3)*y(4)**3 &
1541 & + 831*y(4)**4))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) + y(3) &
1544 y = s_cb(i - 1:i + 2) - s_cb(i - 2:i + 1)
1546 & 0) = (4*y(3)**2*(5*y(1)**2*y(2)**2 + 5*y(1)**2*y(2)*y(3) + 50*y(1)**2*y(3)**2 &
1547 & + 10*y(1)*y(2)**3 + 15*y(1)*y(2)**2*y(3) + 205*y(1)*y(2)*y(3)**2 + 100*y(1)*y(3) &
1548 & **3 + 5*y(2)**4 + 10*y(2)**3*y(3) + 205*y(2)**2*y(3)**2 + 200*y(2)*y(3)**3 &
1549 & + 831*y(3)**4))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) &
1552 & 1) = (4*y(3)**2*(5*y(1)**3*y(2)*y(3) + 10*y(1)**3*y(2)*y(4) - 95*y(1)**3*y(3)**2 &
1553 & + 5*y(1)**3*y(3)*y(4) + 20*y(1)**2*y(2)**2*y(3) + 40*y(1)**2*y(2)**2*y(4) &
1554 & - 465*y(1)**2*y(2)*y(3)**2 + 55*y(1)**2*y(2)*y(3)*y(4) + 10*y(1)**2*y(2)*y(4)**2 &
1555 & - 285*y(1)**2*y(3)**3 + 20*y(1)**2*y(3)**2*y(4) + 5*y(1)**2*y(3)*y(4)**2 &
1556 & + 30*y(1)*y(2)**3*y(3) + 60*y(1)*y(2)**3*y(4) - 825*y(1)*y(2)**2*y(3)**2 &
1557 & + 135*y(1)*y(2)**2*y(3)*y(4) + 30*y(1)*y(2)**2*y(4)**2 - 1040*y(1)*y(2)*y(3)**3 &
1558 & + 100*y(1)*y(2)*y(3)**2*y(4) + 35*y(1)*y(2)*y(3)*y(4)**2 - 1847*y(1)*y(3)**4 &
1559 & + 125*y(1)*y(3)**3*y(4) + 110*y(1)*y(3)**2*y(4)**2 + 15*y(2)**4*y(3) + 30*y(2) &
1560 & **4*y(4) - 550*y(2)**3*y(3)**2 + 90*y(2)**3*y(3)*y(4) + 20*y(2)**3*y(4)**2 &
1561 & - 1040*y(2)**2*y(3)**3 + 100*y(2)**2*y(3)**2*y(4) + 35*y(2)**2*y(3)*y(4)**2 &
1562 & - 3694*y(2)*y(3)**4 + 250*y(2)*y(3)**3*y(4) + 220*y(2)*y(3)**2*y(4)**2 &
1563 & - 3219*y(3)**5 - 1452*y(3)**4*y(4) + 105*y(3)**3*y(4)**2))/(5*(y(2) + y(3))*(y(3) &
1564 & + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4)) &
1567 & 2) = -(4*y(3)**2*(5*y(2)**3*y(3) - 95*y(2)*y(3)**3 - 190*y(2)**2*y(3)**2 &
1568 & + 10*y(2)**3*y(4) + 100*y(3)**3*y(4) - 1562*y(3)**4 - 95*y(1)*y(2)*y(3)**2 &
1569 & + 5*y(1)*y(2)**2*y(3) + 10*y(1)*y(2)**2*y(4) + 100*y(1)*y(3)**2*y(4) + 205*y(2) &
1570 & *y(3)**2*y(4) + 15*y(2)**2*y(3)*y(4) + 10*y(1)*y(2)*y(3)*y(4)))/(5*(y(1) + y(2)) &
1571 & *(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
1574 & 3) = (4*y(3)**2*(50*y(1)**4*y(3)**2 + 5*y(1)**4*y(3)*y(4) + 5*y(1)**4*y(4)**2 &
1575 & + 300*y(1)**3*y(2)*y(3)**2 + 30*y(1)**3*y(2)*y(3)*y(4) + 30*y(1)**3*y(2)*y(4)**2 &
1576 & + 200*y(1)**3*y(3)**3 + 25*y(1)**3*y(3)**2*y(4) + 35*y(1)**3*y(3)*y(4)**2 &
1577 & + 10*y(1)**3*y(4)**3 + 750*y(1)**2*y(2)**2*y(3)**2 + 75*y(1)**2*y(2)**2*y(3)*y(4) &
1578 & + 75*y(1)**2*y(2)**2*y(4)**2 + 1000*y(1)**2*y(2)*y(3)**3 + 125*y(1)**2*y(2)*y(3) &
1579 & **2*y(4) + 175*y(1)**2*y(2)*y(3)*y(4)**2 + 50*y(1)**2*y(2)*y(4)**3 + 1081*y(1) &
1580 & **2*y(3)**4 - 50*y(1)**2*y(3)**3*y(4) - 10*y(1)**2*y(3)**2*y(4)**2 + 45*y(1) &
1581 & **2*y(3)*y(4)**3 + 5*y(1)**2*y(4)**4 + 900*y(1)*y(2)**3*y(3)**2 + 90*y(1)*y(2) &
1582 & **3*y(3)*y(4) + 90*y(1)*y(2)**3*y(4)**2 + 1800*y(1)*y(2)**2*y(3)**3 + 225*y(1) &
1583 & *y(2)**2*y(3)**2*y(4) + 315*y(1)*y(2)**2*y(3)*y(4)**2 + 90*y(1)*y(2)**2*y(4)**3 &
1584 & + 4224*y(1)*y(2)*y(3)**4 - 120*y(1)*y(2)*y(3)**3*y(4) + 25*y(1)*y(2)*y(3)**2*y(4) &
1585 & **2 + 165*y(1)*y(2)*y(3)*y(4)**3 + 20*y(1)*y(2)*y(4)**4 + 3324*y(1)*y(3)**5 &
1586 & + 1407*y(1)*y(3)**4*y(4) - 100*y(1)*y(3)**3*y(4)**2 + 70*y(1)*y(3)**2*y(4)**3 &
1587 & + 15*y(1)*y(3)*y(4)**4 + 450*y(2)**4*y(3)**2 + 45*y(2)**4*y(3)*y(4) + 45*y(2) &
1588 & **4*y(4)**2 + 1200*y(2)**3*y(3)**3 + 150*y(2)**3*y(3)**2*y(4) + 210*y(2)**3*y(3) &
1589 & *y(4)**2 + 60*y(2)**3*y(4)**3 + 4224*y(2)**2*y(3)**4 - 120*y(2)**2*y(3)**3*y(4) &
1590 & + 25*y(2)**2*y(3)**2*y(4)**2 + 165*y(2)**2*y(3)*y(4)**3 + 20*y(2)**2*y(4)**4 &
1591 & + 6648*y(2)*y(3)**5 + 2814*y(2)*y(3)**4*y(4) - 200*y(2)*y(3)**3*y(4)**2 &
1592 & + 140*y(2)*y(3)**2*y(4)**3 + 30*y(2)*y(3)*y(4)**4 + 3174*y(3)**6 + 3039*y(3) &
1593 & **5*y(4) + 771*y(3)**4*y(4)**2 + 135*y(3)**3*y(4)**3 + 60*y(3)**2*y(4)**4)) &
1594 & /(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
1595 & + y(3) + y(4))**2)
1597 & 4) = -(4*y(3)**2*(100*y(1)**2*y(2)*y(3)**2 + 10*y(1)**2*y(2)*y(3)*y(4) + 10*y(1) &
1598 & **2*y(2)*y(4)**2 - 95*y(1)**2*y(3)**2*y(4) + 5*y(1)**2*y(3)*y(4)**2 + 300*y(1) &
1599 & *y(2)**2*y(3)**2 + 30*y(1)*y(2)**2*y(3)*y(4) + 30*y(1)*y(2)**2*y(4)**2 + 200*y(1) &
1600 & *y(2)*y(3)**3 - 260*y(1)*y(2)*y(3)**2*y(4) + 50*y(1)*y(2)*y(3)*y(4)**2 + 10*y(1) &
1601 & *y(2)*y(4)**3 + 1562*y(1)*y(3)**4 - 190*y(1)*y(3)**3*y(4) + 15*y(1)*y(3)**2*y(4) &
1602 & **2 + 5*y(1)*y(3)*y(4)**3 + 300*y(2)**3*y(3)**2 + 30*y(2)**3*y(3)*y(4) + 30*y(2) &
1603 & **3*y(4)**2 + 400*y(2)**2*y(3)**3 - 235*y(2)**2*y(3)**2*y(4) + 85*y(2)**2*y(3) &
1604 & *y(4)**2 + 20*y(2)**2*y(4)**3 + 3224*y(2)*y(3)**4 - 460*y(2)*y(3)**3*y(4) &
1605 & - 35*y(2)*y(3)**2*y(4)**2 + 25*y(2)*y(3)*y(4)**3 + 3124*y(3)**5 + 1467*y(3) &
1606 & **4*y(4) + 110*y(3)**3*y(4)**2 + 105*y(3)**2*y(4)**3))/(5*(y(1) + y(2))*(y(2) &
1607 & + y(3))*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)) &
1610 & 5) = (4*y(3)**2*(50*y(2)**2*y(3)**2 + 5*y(2)**2*y(3)*y(4) + 5*y(2)**2*y(4)**2 &
1611 & - 95*y(2)*y(3)**2*y(4) + 5*y(2)*y(3)*y(4)**2 + 781*y(3)**4 + 50*y(3)**2*y(4)**2)) &
1612 & /(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) + y(3) + y(4))**2)
1614 y = s_cb(i:i + 3) - s_cb(i - 1:i + 2)
1616 & 0) = (4*y(2)**2*(50*y(1)**2*y(2)**2 + 5*y(1)**2*y(2)*y(3) + 5*y(1)**2*y(3)**2 &
1617 & - 95*y(1)*y(2)**2*y(3) + 5*y(1)*y(2)*y(3)**2 + 781*y(2)**4 + 50*y(2)**2*y(3)**2)) &
1618 & /(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
1620 & 1) = -(4*y(2)**2*(105*y(1)**3*y(2)**2 + 25*y(1)**3*y(2)*y(3) + 5*y(1)**3*y(2) &
1621 & *y(4) + 20*y(1)**3*y(3)**2 + 10*y(1)**3*y(3)*y(4) + 110*y(1)**2*y(2)**3 - 35*y(1) &
1622 & **2*y(2)**2*y(3) + 15*y(1)**2*y(2)**2*y(4) + 85*y(1)**2*y(2)*y(3)**2 + 50*y(1) &
1623 & **2*y(2)*y(3)*y(4) + 5*y(1)**2*y(2)*y(4)**2 + 30*y(1)**2*y(3)**3 + 30*y(1) &
1624 & **2*y(3)**2*y(4) + 10*y(1)**2*y(3)*y(4)**2 + 1467*y(1)*y(2)**4 - 460*y(1)*y(2) &
1625 & **3*y(3) - 190*y(1)*y(2)**3*y(4) - 235*y(1)*y(2)**2*y(3)**2 - 260*y(1)*y(2) &
1626 & **2*y(3)*y(4) - 95*y(1)*y(2)**2*y(4)**2 + 30*y(1)*y(2)*y(3)**3 + 30*y(1)*y(2) &
1627 & *y(3)**2*y(4) + 10*y(1)*y(2)*y(3)*y(4)**2 + 3124*y(2)**5 + 3224*y(2)**4*y(3) &
1628 & + 1562*y(2)**4*y(4) + 400*y(2)**3*y(3)**2 + 200*y(2)**3*y(3)*y(4) + 300*y(2) &
1629 & **2*y(3)**3 + 300*y(2)**2*y(3)**2*y(4) + 100*y(2)**2*y(3)*y(4)**2))/(5*(y(2) &
1630 & + y(3))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
1631 & + y(3) + y(4))**2)
1633 & 2) = -(4*y(2)**2*(100*y(1)*y(2)**3 - 190*y(2)**2*y(3)**2 + 10*y(1)*y(3)**3 &
1634 & + 5*y(2)*y(3)**3 - 95*y(2)**3*y(3) - 1562*y(2)**4 + 15*y(1)*y(2)*y(3)**2 &
1635 & + 205*y(1)*y(2)**2*y(3) + 100*y(1)*y(2)**2*y(4) + 10*y(1)*y(3)**2*y(4) + 5*y(2) &
1636 & *y(3)**2*y(4) - 95*y(2)**2*y(3)*y(4) + 10*y(1)*y(2)*y(3)*y(4)))/(5*(y(1) + y(2)) &
1637 & *(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
1640 & 3) = (4*y(2)**2*(60*y(1)**4*y(2)**2 + 30*y(1)**4*y(2)*y(3) + 15*y(1)**4*y(2)*y(4) &
1641 & + 20*y(1)**4*y(3)**2 + 20*y(1)**4*y(3)*y(4) + 5*y(1)**4*y(4)**2 + 135*y(1) &
1642 & **3*y(2)**3 + 140*y(1)**3*y(2)**2*y(3) + 70*y(1)**3*y(2)**2*y(4) + 165*y(1) &
1643 & **3*y(2)*y(3)**2 + 165*y(1)**3*y(2)*y(3)*y(4) + 45*y(1)**3*y(2)*y(4)**2 + 60*y(1) &
1644 & **3*y(3)**3 + 90*y(1)**3*y(3)**2*y(4) + 50*y(1)**3*y(3)*y(4)**2 + 10*y(1)**3*y(4) &
1645 & **3 + 771*y(1)**2*y(2)**4 - 200*y(1)**2*y(2)**3*y(3) - 100*y(1)**2*y(2)**3*y(4) &
1646 & + 25*y(1)**2*y(2)**2*y(3)**2 + 25*y(1)**2*y(2)**2*y(3)*y(4) - 10*y(1)**2*y(2) &
1647 & **2*y(4)**2 + 210*y(1)**2*y(2)*y(3)**3 + 315*y(1)**2*y(2)*y(3)**2*y(4) + 175*y(1) &
1648 & **2*y(2)*y(3)*y(4)**2 + 35*y(1)**2*y(2)*y(4)**3 + 45*y(1)**2*y(3)**4 + 90*y(1) &
1649 & **2*y(3)**3*y(4) + 75*y(1)**2*y(3)**2*y(4)**2 + 30*y(1)**2*y(3)*y(4)**3 + 5*y(1) &
1650 & **2*y(4)**4 + 3039*y(1)*y(2)**5 + 2814*y(1)*y(2)**4*y(3) + 1407*y(1)*y(2)**4*y(4) &
1651 & - 120*y(1)*y(2)**3*y(3)**2 - 120*y(1)*y(2)**3*y(3)*y(4) - 50*y(1)*y(2)**3*y(4) &
1652 & **2 + 150*y(1)*y(2)**2*y(3)**3 + 225*y(1)*y(2)**2*y(3)**2*y(4) + 125*y(1)*y(2) &
1653 & **2*y(3)*y(4)**2 + 25*y(1)*y(2)**2*y(4)**3 + 45*y(1)*y(2)*y(3)**4 + 90*y(1)*y(2) &
1654 & *y(3)**3*y(4) + 75*y(1)*y(2)*y(3)**2*y(4)**2 + 30*y(1)*y(2)*y(3)*y(4)**3 + 5*y(1) &
1655 & *y(2)*y(4)**4 + 3174*y(2)**6 + 6648*y(2)**5*y(3) + 3324*y(2)**5*y(4) + 4224*y(2) &
1656 & **4*y(3)**2 + 4224*y(2)**4*y(3)*y(4) + 1081*y(2)**4*y(4)**2 + 1200*y(2)**3*y(3) &
1657 & **3 + 1800*y(2)**3*y(3)**2*y(4) + 1000*y(2)**3*y(3)*y(4)**2 + 200*y(2)**3*y(4) &
1658 & **3 + 450*y(2)**2*y(3)**4 + 900*y(2)**2*y(3)**3*y(4) + 750*y(2)**2*y(3)**2*y(4) &
1659 & **2 + 300*y(2)**2*y(3)*y(4)**3 + 50*y(2)**2*y(4)**4))/(5*(y(2) + y(3))**2*(y(1) &
1660 & + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
1662 & 4) = (4*y(2)**2*(105*y(1)**2*y(2)**3 + 220*y(1)**2*y(2)**2*y(3) + 110*y(1) &
1663 & **2*y(2)**2*y(4) + 35*y(1)**2*y(2)*y(3)**2 + 35*y(1)**2*y(2)*y(3)*y(4) + 5*y(1) &
1664 & **2*y(2)*y(4)**2 + 20*y(1)**2*y(3)**3 + 30*y(1)**2*y(3)**2*y(4) + 10*y(1)**2*y(3) &
1665 & *y(4)**2 - 1452*y(1)*y(2)**4 + 250*y(1)*y(2)**3*y(3) + 125*y(1)*y(2)**3*y(4) &
1666 & + 100*y(1)*y(2)**2*y(3)**2 + 100*y(1)*y(2)**2*y(3)*y(4) + 20*y(1)*y(2)**2*y(4) &
1667 & **2 + 90*y(1)*y(2)*y(3)**3 + 135*y(1)*y(2)*y(3)**2*y(4) + 55*y(1)*y(2)*y(3)*y(4) &
1668 & **2 + 5*y(1)*y(2)*y(4)**3 + 30*y(1)*y(3)**4 + 60*y(1)*y(3)**3*y(4) + 40*y(1)*y(3) &
1669 & **2*y(4)**2 + 10*y(1)*y(3)*y(4)**3 - 3219*y(2)**5 - 3694*y(2)**4*y(3) - 1847*y(2) &
1670 & **4*y(4) - 1040*y(2)**3*y(3)**2 - 1040*y(2)**3*y(3)*y(4) - 285*y(2)**3*y(4)**2 &
1671 & - 550*y(2)**2*y(3)**3 - 825*y(2)**2*y(3)**2*y(4) - 465*y(2)**2*y(3)*y(4)**2 &
1672 & - 95*y(2)**2*y(4)**3 + 15*y(2)*y(3)**4 + 30*y(2)*y(3)**3*y(4) + 20*y(2)*y(3) &
1673 & **2*y(4)**2 + 5*y(2)*y(3)*y(4)**3))/(5*(y(1) + y(2))*(y(2) + y(3))*(y(1) + y(2) &
1674 & + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
1676 & 5) = (4*y(2)**2*(831*y(2)**4 + 200*y(2)**3*y(3) + 100*y(2)**3*y(4) + 205*y(2) &
1677 & **2*y(3)**2 + 205*y(2)**2*y(3)*y(4) + 50*y(2)**2*y(4)**2 + 10*y(2)*y(3)**3 &
1678 & + 15*y(2)*y(3)**2*y(4) + 5*y(2)*y(3)*y(4)**2 + 5*y(3)**4 + 10*y(3)**3*y(4) &
1679 & + 5*y(3)**2*y(4)**2))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) &
1680 & + y(3) + y(4))**2)
1682 y = s_cb(i + 1:i + 4) - s_cb(i:i + 3)
1684 & 0) = (4*y(1)**2*(831*y(1)**4 + 200*y(1)**3*y(2) + 100*y(1)**3*y(3) + 205*y(1) &
1685 & **2*y(2)**2 + 205*y(1)**2*y(2)*y(3) + 50*y(1)**2*y(3)**2 + 10*y(1)*y(2)**3 &
1686 & + 15*y(1)*y(2)**2*y(3) + 5*y(1)*y(2)*y(3)**2 + 5*y(2)**4 + 10*y(2)**3*y(3) &
1687 & + 5*y(2)**2*y(3)**2))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
1688 & + y(3) + y(4))**2)
1690 & 1) = -(4*y(1)**2*(1662*y(1)**5 + 3824*y(1)**4*y(2) + 3624*y(1)**4*y(3) &
1691 & + 1762*y(1)**4*y(4) + 1515*y(1)**3*y(2)**2 + 2115*y(1)**3*y(2)*y(3) + 805*y(1) &
1692 & **3*y(2)*y(4) + 700*y(1)**3*y(3)**2 + 500*y(1)**3*y(3)*y(4) + 100*y(1)**3*y(4) &
1693 & **2 + 1060*y(1)**2*y(2)**3 + 2205*y(1)**2*y(2)**2*y(3) + 835*y(1)**2*y(2)**2*y(4) &
1694 & + 1445*y(1)**2*y(2)*y(3)**2 + 1030*y(1)**2*y(2)*y(3)*y(4) + 205*y(1)**2*y(2)*y(4) &
1695 & **2 + 300*y(1)**2*y(3)**3 + 300*y(1)**2*y(3)**2*y(4) + 100*y(1)**2*y(3)*y(4)**2 &
1696 & + 75*y(1)*y(2)**4 + 180*y(1)*y(2)**3*y(3) + 60*y(1)*y(2)**3*y(4) + 135*y(1)*y(2) &
1697 & **2*y(3)**2 + 90*y(1)*y(2)**2*y(3)*y(4) + 15*y(1)*y(2)**2*y(4)**2 + 30*y(1)*y(2) &
1698 & *y(3)**3 + 30*y(1)*y(2)*y(3)**2*y(4) + 10*y(1)*y(2)*y(3)*y(4)**2 + 30*y(2)**5 &
1699 & + 90*y(2)**4*y(3) + 30*y(2)**4*y(4) + 90*y(2)**3*y(3)**2 + 60*y(2)**3*y(3)*y(4) &
1700 & + 10*y(2)**3*y(4)**2 + 30*y(2)**2*y(3)**3 + 30*y(2)**2*y(3)**2*y(4) + 10*y(2) &
1701 & **2*y(3)*y(4)**2))/(5*(y(2) + y(3))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) &
1702 & + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
1704 & 2) = (4*y(1)**2*(1767*y(1)**4 + 725*y(1)**3*y(2) + 415*y(1)**3*y(3) + 105*y(4) &
1705 & *y(1)**3 + 665*y(1)**2*y(2)**2 + 775*y(1)**2*y(2)*y(3) + 220*y(4)*y(1)**2*y(2) &
1706 & + 215*y(1)**2*y(3)**2 + 110*y(4)*y(1)**2*y(3) + 75*y(1)*y(2)**3 + 130*y(1)*y(2) &
1707 & **2*y(3) + 35*y(4)*y(1)*y(2)**2 + 60*y(1)*y(2)*y(3)**2 + 35*y(4)*y(1)*y(2)*y(3) &
1708 & + 5*y(1)*y(3)**3 + 5*y(4)*y(1)*y(3)**2 + 30*y(2)**4 + 70*y(2)**3*y(3) + 20*y(4) &
1709 & *y(2)**3 + 50*y(2)**2*y(3)**2 + 30*y(4)*y(2)**2*y(3) + 10*y(2)*y(3)**3 + 10*y(4) &
1710 & *y(2)*y(3)**2))/(5*(y(1) + y(2))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) &
1711 & + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
1713 & 3) = (4*y(1)**2*(831*y(1)**6 + 3624*y(1)**5*y(2) + 3524*y(1)**5*y(3) + 1762*y(1) &
1714 & **5*y(4) + 4884*y(1)**4*y(2)**2 + 9058*y(1)**4*y(2)*y(3) + 4529*y(1)**4*y(2)*y(4) &
1715 & + 4224*y(1)**4*y(3)**2 + 4224*y(1)**4*y(3)*y(4) + 1081*y(1)**4*y(4)**2 &
1716 & + 2565*y(1)**3*y(2)**3 + 6120*y(1)**3*y(2)**2*y(3) + 3060*y(1)**3*y(2)**2*y(4) &
1717 & + 4755*y(1)**3*y(2)*y(3)**2 + 4755*y(1)**3*y(2)*y(3)*y(4) + 1315*y(1)**3*y(2) &
1718 & *y(4)**2 + 1200*y(1)**3*y(3)**3 + 1800*y(1)**3*y(3)**2*y(4) + 1000*y(1)**3*y(3) &
1719 & *y(4)**2 + 200*y(1)**3*y(4)**3 + 1395*y(1)**2*y(2)**4 + 4380*y(1)**2*y(2)**3*y(3) &
1720 & + 2190*y(1)**2*y(2)**3*y(4) + 5025*y(1)**2*y(2)**2*y(3)**2 + 5025*y(1)**2*y(2) &
1721 & **2*y(3)*y(4) + 1390*y(1)**2*y(2)**2*y(4)**2 + 2490*y(1)**2*y(2)*y(3)**3 &
1722 & + 3735*y(1)**2*y(2)*y(3)**2*y(4) + 2075*y(1)**2*y(2)*y(3)*y(4)**2 + 415*y(1) &
1723 & **2*y(2)*y(4)**3 + 450*y(1)**2*y(3)**4 + 900*y(1)**2*y(3)**3*y(4) + 750*y(1) &
1724 & **2*y(3)**2*y(4)**2 + 300*y(1)**2*y(3)*y(4)**3 + 50*y(1)**2*y(4)**4 + 135*y(1) &
1725 & *y(2)**5 + 450*y(1)*y(2)**4*y(3) + 225*y(1)*y(2)**4*y(4) + 540*y(1)*y(2)**3*y(3) &
1726 & **2 + 540*y(1)*y(2)**3*y(3)*y(4) + 150*y(1)*y(2)**3*y(4)**2 + 270*y(1)*y(2) &
1727 & **2*y(3)**3 + 405*y(1)*y(2)**2*y(3)**2*y(4) + 225*y(1)*y(2)**2*y(3)*y(4)**2 &
1728 & + 45*y(1)*y(2)**2*y(4)**3 + 45*y(1)*y(2)*y(3)**4 + 90*y(1)*y(2)*y(3)**3*y(4) &
1729 & + 75*y(1)*y(2)*y(3)**2*y(4)**2 + 30*y(1)*y(2)*y(3)*y(4)**3 + 5*y(1)*y(2)*y(4)**4 &
1730 & + 45*y(2)**6 + 180*y(2)**5*y(3) + 90*y(2)**5*y(4) + 270*y(2)**4*y(3)**2 &
1731 & + 270*y(2)**4*y(3)*y(4) + 75*y(2)**4*y(4)**2 + 180*y(2)**3*y(3)**3 + 270*y(2) &
1732 & **3*y(3)**2*y(4) + 150*y(2)**3*y(3)*y(4)**2 + 30*y(2)**3*y(4)**3 + 45*y(2) &
1733 & **2*y(3)**4 + 90*y(2)**2*y(3)**3*y(4) + 75*y(2)**2*y(3)**2*y(4)**2 + 30*y(2) &
1734 & **2*y(3)*y(4)**3 + 5*y(2)**2*y(4)**4))/(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3)) &
1735 & **2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
1737 & 4) = -(4*y(1)**2*(1767*y(1)**5 + 4464*y(1)**4*y(2) + 4154*y(1)**4*y(3) &
1738 & + 2077*y(1)**4*y(4) + 2655*y(1)**3*y(2)**2 + 4010*y(1)**3*y(2)*y(3) + 2005*y(1) &
1739 & **3*y(2)*y(4) + 1460*y(1)**3*y(3)**2 + 1460*y(1)**3*y(3)*y(4) + 415*y(1)**3*y(4) &
1740 & **2 + 1800*y(1)**2*y(2)**3 + 4000*y(1)**2*y(2)**2*y(3) + 2000*y(1)**2*y(2) &
1741 & **2*y(4) + 2850*y(1)**2*y(2)*y(3)**2 + 2850*y(1)**2*y(2)*y(3)*y(4) + 790*y(1) &
1742 & **2*y(2)*y(4)**2 + 650*y(1)**2*y(3)**3 + 975*y(1)**2*y(3)**2*y(4) + 535*y(1) &
1743 & **2*y(3)*y(4)**2 + 105*y(1)**2*y(4)**3 + 270*y(1)*y(2)**4 + 720*y(1)*y(2)**3*y(3) &
1744 & + 360*y(1)*y(2)**3*y(4) + 645*y(1)*y(2)**2*y(3)**2 + 645*y(1)*y(2)**2*y(3)*y(4) &
1745 & + 165*y(1)*y(2)**2*y(4)**2 + 210*y(1)*y(2)*y(3)**3 + 315*y(1)*y(2)*y(3)**2*y(4) &
1746 & + 155*y(1)*y(2)*y(3)*y(4)**2 + 25*y(1)*y(2)*y(4)**3 + 15*y(1)*y(3)**4 + 30*y(1) &
1747 & *y(3)**3*y(4) + 20*y(1)*y(3)**2*y(4)**2 + 5*y(1)*y(3)*y(4)**3 + 90*y(2)**5 &
1748 & + 300*y(2)**4*y(3) + 150*y(2)**4*y(4) + 360*y(2)**3*y(3)**2 + 360*y(2)**3*y(3) &
1749 & *y(4) + 90*y(2)**3*y(4)**2 + 180*y(2)**2*y(3)**3 + 270*y(2)**2*y(3)**2*y(4) &
1750 & + 130*y(2)**2*y(3)*y(4)**2 + 20*y(2)**2*y(4)**3 + 30*y(2)*y(3)**4 + 60*y(2)*y(3) &
1751 & **3*y(4) + 40*y(2)*y(3)**2*y(4)**2 + 10*y(2)*y(3)*y(4)**3))/(5*(y(1) + y(2)) &
1752 & *(y(2) + y(3))*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
1755 & 5) = (4*y(1)**2*(996*y(1)**4 + 675*y(1)**3*y(2) + 450*y(1)**3*y(3) + 225*y(1) &
1756 & **3*y(4) + 600*y(1)**2*y(2)**2 + 800*y(1)**2*y(2)*y(3) + 400*y(1)**2*y(2)*y(4) &
1757 & + 260*y(1)**2*y(3)**2 + 260*y(1)**2*y(3)*y(4) + 60*y(1)**2*y(4)**2 + 135*y(1) &
1758 & *y(2)**3 + 270*y(1)*y(2)**2*y(3) + 135*y(1)*y(2)**2*y(4) + 165*y(1)*y(2)*y(3)**2 &
1759 & + 165*y(1)*y(2)*y(3)*y(4) + 30*y(1)*y(2)*y(4)**2 + 30*y(1)*y(3)**3 + 45*y(1)*y(3) &
1760 & **2*y(4) + 15*y(1)*y(3)*y(4)**2 + 45*y(2)**4 + 120*y(2)**3*y(3) + 60*y(2)**3*y(4) &
1761 & + 110*y(2)**2*y(3)**2 + 110*y(2)**2*y(3)*y(4) + 20*y(2)**2*y(4)**2 + 40*y(2)*y(3) &
1762 & **3 + 60*y(2)*y(3)**2*y(4) + 20*y(2)*y(3)*y(4)**2 + 5*y(3)**4 + 10*y(3)**3*y(4) &
1763 & + 5*y(3)**2*y(4)**2))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) &
1764 & + y(3) + y(4))**2)
1782# 192 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
1784 if (weno_dir == 2)
then
1785 if (weno_order == 3)
then
1786 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
1788 poly_coef_cbr_y(i + 1, 0, 0) = (s_cb(i) - s_cb(i + 1))/(s_cb(i) - s_cb(i + 2))
1789 poly_coef_cbr_y(i + 1, 1, 0) = (s_cb(i) - s_cb(i + 1))/(s_cb(i - 1) - s_cb(i + 1))
1795 d_cbr_y(0, i + 1) = (s_cb(i - 1) - s_cb(i + 1))/(s_cb(i - 1) - s_cb(i + 2))
1796 d_cbl_y(0, i + 1) = (s_cb(i - 1) - s_cb(i))/(s_cb(i - 1) - s_cb(i + 2))
1802 beta_coef_y(i + 1, 0, 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp/(s_cb(i) - s_cb(i + 2))**2._wp
1803 beta_coef_y(i + 1, 1, 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp/(s_cb(i - 1) - s_cb(i + 1))**2._wp
1808 if (null_weights)
then
1809 if (bc_s%beg == bc_riemann_extrap)
then
1814 if (bc_s%end == bc_riemann_extrap)
then
1822 else if (weno_order == 5)
then
1823 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
1826 & 0) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i + 2)))/((s_cb(i) - s_cb(i &
1827 & + 3))*(s_cb(i + 3) - s_cb(i + 1)))
1829 & 0) = ((s_cb(i - 1) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i)))/((s_cb(i - 1) &
1830 & - s_cb(i + 2))*(s_cb(i + 2) - s_cb(i)))
1832 & 1) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i + 2)))/((s_cb(i - 1) &
1833 & - s_cb(i + 1))*(s_cb(i - 1) - s_cb(i + 2)))
1835 & 1) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))/((s_cb(i - 2) &
1836 & - s_cb(i))*(s_cb(i - 2) - s_cb(i + 1)))
1838 & 0) = ((s_cb(i + 1) - s_cb(i))*(s_cb(i) - s_cb(i + 2)))/((s_cb(i) - s_cb(i + 3)) &
1839 & *(s_cb(i + 3) - s_cb(i + 1)))
1841 & 0) = ((s_cb(i) - s_cb(i - 1))*(s_cb(i) - s_cb(i + 1)))/((s_cb(i - 1) - s_cb(i &
1842 & + 2))*(s_cb(i) - s_cb(i + 2)))
1844 & 1) = ((s_cb(i + 1) - s_cb(i))*(s_cb(i) - s_cb(i + 2)))/((s_cb(i - 1) - s_cb(i &
1845 & + 1))*(s_cb(i - 1) - s_cb(i + 2)))
1847 & 1) = ((s_cb(i - 1) - s_cb(i))*(s_cb(i) - s_cb(i + 1)))/((s_cb(i - 2) - s_cb(i)) &
1848 & *(s_cb(i - 2) - s_cb(i + 1)))
1851 & 1) = ((s_cb(i) - s_cb(i + 2)) + (s_cb(i + 1) - s_cb(i + 3)))/((s_cb(i) - s_cb(i &
1852 & + 2))*(s_cb(i) - s_cb(i + 3)))*((s_cb(i) - s_cb(i + 1)))
1854 & 0) = ((s_cb(i - 2) - s_cb(i + 1)) + (s_cb(i - 1) - s_cb(i + 1)))/((s_cb(i - 1) &
1855 & - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 2)))*((s_cb(i + 1) - s_cb(i)))
1857 & 1) = ((s_cb(i) - s_cb(i + 2)) + (s_cb(i) - s_cb(i + 3)))/((s_cb(i) - s_cb(i + 2)) &
1858 & *(s_cb(i) - s_cb(i + 3)))*((s_cb(i + 1) - s_cb(i)))
1860 & 0) = ((s_cb(i - 2) - s_cb(i)) + (s_cb(i - 1) - s_cb(i + 1)))/((s_cb(i - 2) &
1861 & - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))*((s_cb(i) - s_cb(i + 1)))
1865 & i + 1) = ((s_cb(i - 2) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))/((s_cb(i - 2) &
1866 & - s_cb(i + 3))*(s_cb(i + 3) - s_cb(i - 1)))
1868 & i + 1) = ((s_cb(i + 1) - s_cb(i + 2))*(s_cb(i + 1) - s_cb(i + 3)))/((s_cb(i - 2) &
1869 & - s_cb(i + 2))*(s_cb(i - 2) - s_cb(i + 3)))
1871 & i + 1) = ((s_cb(i - 2) - s_cb(i))*(s_cb(i) - s_cb(i - 1)))/((s_cb(i - 2) - s_cb(i + 3)) &
1872 & *(s_cb(i + 3) - s_cb(i - 1)))
1874 & i + 1) = ((s_cb(i) - s_cb(i + 2))*(s_cb(i) - s_cb(i + 3)))/((s_cb(i - 2) - s_cb(i + 2)) &
1875 & *(s_cb(i - 2) - s_cb(i + 3)))
1882 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1883 & + (s_cb(i + 1) - s_cb(i))*(s_cb(i + 2) - s_cb(i + 1)) + (s_cb(i + 2) - s_cb(i + 1)) &
1884 & **2._wp)/((s_cb(i) - s_cb(i + 3))**2._wp*(s_cb(i + 1) - s_cb(i + 3))**2._wp)
1887 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(19._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1888 & - (s_cb(i + 1) - s_cb(i))*(s_cb(i + 3) - s_cb(i + 1)) + 2._wp*(s_cb(i + 2) - s_cb(i)) &
1889 & *((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1))))/((s_cb(i) - s_cb(i + 2)) &
1890 & *(s_cb(i) - s_cb(i + 3))**2._wp*(s_cb(i + 3) - s_cb(i + 1)))
1893 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1894 & + (s_cb(i + 1) - s_cb(i))*((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1))) &
1895 & + ((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1)))**2._wp)/((s_cb(i) - s_cb(i &
1896 & + 2))**2._wp*(s_cb(i) - s_cb(i + 3))**2._wp)
1899 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1900 & + (s_cb(i) - s_cb(i - 1))**2._wp + (s_cb(i) - s_cb(i - 1))*(s_cb(i + 1) - s_cb(i))) &
1901 & /((s_cb(i - 1) - s_cb(i + 2))**2._wp*(s_cb(i) - s_cb(i + 2))**2._wp)
1904 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*((s_cb(i) - s_cb(i + 1))*((s_cb(i) &
1905 & - s_cb(i - 1)) + 20._wp*(s_cb(i + 1) - s_cb(i))) + (2._wp*(s_cb(i) - s_cb(i - 1)) &
1906 & + (s_cb(i + 1) - s_cb(i)))*(s_cb(i + 2) - s_cb(i)))/((s_cb(i + 1) - s_cb(i - 1)) &
1907 & *(s_cb(i - 1) - s_cb(i + 2))**2._wp*(s_cb(i + 2) - s_cb(i)))
1910 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1911 & + (s_cb(i + 1) - s_cb(i))*(s_cb(i + 2) - s_cb(i + 1)) + (s_cb(i + 2) - s_cb(i + 1)) &
1912 & **2._wp)/((s_cb(i - 1) - s_cb(i + 1))**2._wp*(s_cb(i - 1) - s_cb(i + 2))**2._wp)
1915 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(12._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1916 & + ((s_cb(i) - s_cb(i - 2)) + (s_cb(i) - s_cb(i - 1)))**2._wp + 3._wp*((s_cb(i) &
1917 & - s_cb(i - 2)) + (s_cb(i) - s_cb(i - 1)))*(s_cb(i + 1) - s_cb(i)))/((s_cb(i - 2) &
1918 & - s_cb(i + 1))**2._wp*(s_cb(i - 1) - s_cb(i + 1))**2._wp)
1921 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(19._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1922 & + ((s_cb(i) - s_cb(i - 2))*(s_cb(i) - s_cb(i + 1))) + 2._wp*(s_cb(i + 1) - s_cb(i &
1923 & - 1))*((s_cb(i) - s_cb(i - 2)) + (s_cb(i + 1) - s_cb(i - 1))))/((s_cb(i - 2) &
1924 & - s_cb(i))*(s_cb(i - 2) - s_cb(i + 1))**2._wp*(s_cb(i + 1) - s_cb(i - 1)))
1927 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
1928 & + (s_cb(i) - s_cb(i - 1))**2._wp + (s_cb(i) - s_cb(i - 1))*(s_cb(i + 1) - s_cb(i))) &
1929 & /((s_cb(i - 2) - s_cb(i))**2._wp*(s_cb(i - 2) - s_cb(i + 1))**2._wp)
1934 if (null_weights)
then
1935 if (bc_s%beg == bc_riemann_extrap)
then
1942 if (bc_s%end == bc_riemann_extrap)
then
1944 & s - 1)/sum(
d_cbr_y(:,s - 1))
1946 & s - 1)/sum(
d_cbl_y(:,s - 1))
1952 if (.not. teno)
then
1953 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
1966 w = s_cb(i - 3:i + 4) - s_cb(i)
1968 & i + 1) = ((w(5) - w(6))*(w(5) - w(7))*(w(5) - w(8)))/((w(1) - w(6))*(w(1) - w(7)) &
1971 & i + 1) = ((w(1) - w(5))*(w(5) - w(7))*(w(5) - w(8))*(w(1)*w(2) - w(1)*w(6) - w(1) &
1972 & *w(7) - w(2)*w(6) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) + w(6)*w(7) + w(6)*w(8) + w(7) &
1973 & *w(8) + w(1)**2 + w(2)**2))/((w(1) - w(6))*(w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7)) &
1976 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(5) - w(8))*(w(1)*w(2) + w(1)*w(3) + w(2) &
1977 & *w(3) - w(1)*w(7) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) - w(3)*w(7) - w(3)*w(8) + w(7) &
1978 & *w(8) + w(7)**2 + w(8)**2))/((w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7))*(w(2) - w(8)) &
1981 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(3) - w(5)))/((w(1) - w(8))*(w(2) - w(8)) &
1988 w(1) = s_cb(i + 4) - s_cb(i)
1989 w(2) = s_cb(i + 3) - s_cb(i)
1990 w(3) = s_cb(i + 2) - s_cb(i)
1991 w(4) = s_cb(i + 1) - s_cb(i)
1992 w(5) = s_cb(i) - s_cb(i)
1993 w(6) = s_cb(i - 1) - s_cb(i)
1994 w(7) = s_cb(i - 2) - s_cb(i)
1995 w(8) = s_cb(i - 3) - s_cb(i)
1997 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(3) - w(5)))/((w(1) - w(8))*(w(2) - w(8)) &
2000 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(5) - w(8))*(w(1)*w(2) + w(1)*w(3) + w(2) &
2001 & *w(3) - w(1)*w(7) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) - w(3)*w(7) - w(3)*w(8) + w(7) &
2002 & *w(8) + w(7)**2 + w(8)**2))/((w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7))*(w(2) - w(8)) &
2005 & i + 1) = ((w(1) - w(5))*(w(5) - w(7))*(w(5) - w(8))*(w(1)*w(2) - w(1)*w(6) - w(1) &
2006 & *w(7) - w(2)*w(6) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) + w(6)*w(7) + w(6)*w(8) + w(7) &
2007 & *w(8) + w(1)**2 + w(2)**2))/((w(1) - w(6))*(w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7)) &
2010 & i + 1) = ((w(5) - w(6))*(w(5) - w(7))*(w(5) - w(8)))/((w(1) - w(6))*(w(1) - w(7)) &
2014 y = s_cb(i + 1:i + 4) - s_cb(i:i + 3)
2016 & 0) = (y(1)*y(2)*(y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2017 & + y(2) + y(3) + y(4)))
2019 & 1) = -(y(1)*y(2)*(3*y(2)**2 + 6*y(2)*y(3) + 3*y(2)*y(4) + 2*y(1)*y(2) &
2020 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2021 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2023 & 2) = (y(1)*(y(1)**2 + 3*y(1)*y(2) + 2*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2024 & + 4*y(2)*y(3) + 2*y(4)*y(2) + y(3)**2 + y(4)*y(3)))/((y(1) + y(2))*(y(1) &
2025 & + y(2) + y(3))*(y(1) + y(2) + y(3) + y(4)))
2027 y = s_cb(i:i + 3) - s_cb(i - 1:i + 2)
2029 & 0) = -(y(2)*y(3)*(y(1) + y(2)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2030 & + y(2) + y(3) + y(4)))
2032 & 1) = (y(2)*(y(1) + y(2))*(y(2)**2 + 4*y(2)*y(3) + 2*y(2)*y(4) + y(1)*y(2) &
2033 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2034 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2036 & 2) = (y(2)*y(3)*(y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2037 & + y(2) + y(3) + y(4)))
2039 y = s_cb(i - 1:i + 2) - s_cb(i - 2:i + 1)
2041 & 0) = (y(3)*(y(2) + y(3))*(y(1) + y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) &
2042 & + y(4))*(y(1) + y(2) + y(3) + y(4)))
2044 & 1) = (y(3)*y(4)*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2045 & + 6*y(2)*y(3) + 2*y(4)*y(2) + 3*y(3)**2 + 2*y(4)*y(3)))/((y(2) + y(3))*(y(1) &
2046 & + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2048 & 2) = -(y(3)*y(4)*(y(2) + y(3)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2049 & + y(2) + y(3) + y(4)))
2051 y = s_cb(i - 2:i + 1) - s_cb(i - 3:i)
2053 & 0) = (y(4)*(y(2)**2 + 4*y(2)*y(3) + 4*y(2)*y(4) + y(1)*y(2) + 3*y(3)**2 &
2054 & + 6*y(3)*y(4) + 2*y(1)*y(3) + 3*y(4)**2 + 2*y(1)*y(4)))/((y(3) + y(4))*(y(2) &
2055 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2057 & 1) = -(y(4)*(y(3) + y(4))*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + 2*y(1)*y(4) &
2058 & + 3*y(2)**2 + 6*y(2)*y(3) + 4*y(2)*y(4) + 3*y(3)**2 + 4*y(3)*y(4) + y(4)**2)) &
2059 & /((y(2) + y(3))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2062 & 2) = (y(4)*(y(3) + y(4))*(y(2) + y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) &
2063 & + y(3))*(y(1) + y(2) + y(3) + y(4)))
2066 y(1) = s_cb(i + 1) - s_cb(i)
2067 y(2) = s_cb(i) - s_cb(i - 1)
2068 y(3) = s_cb(i - 1) - s_cb(i - 2)
2069 y(4) = s_cb(i - 2) - s_cb(i - 3)
2071 & 2) = (y(1)*y(2)*(y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2072 & + y(2) + y(3) + y(4)))
2074 & 1) = -(y(1)*y(2)*(3*y(2)**2 + 6*y(2)*y(3) + 3*y(2)*y(4) + 2*y(1)*y(2) &
2075 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2076 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2078 & 0) = (y(1)*(y(1)**2 + 3*y(1)*y(2) + 2*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2079 & + 4*y(2)*y(3) + 2*y(4)*y(2) + y(3)**2 + y(4)*y(3)))/((y(1) + y(2))*(y(1) &
2080 & + y(2) + y(3))*(y(1) + y(2) + y(3) + y(4)))
2083 y(1) = s_cb(i + 2) - s_cb(i + 1)
2084 y(2) = s_cb(i + 1) - s_cb(i)
2085 y(3) = s_cb(i) - s_cb(i - 1)
2086 y(4) = s_cb(i - 1) - s_cb(i - 2)
2088 & 2) = -(y(2)*y(3)*(y(1) + y(2)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2089 & + y(2) + y(3) + y(4)))
2091 & 1) = (y(2)*(y(1) + y(2))*(y(2)**2 + 4*y(2)*y(3) + 2*y(2)*y(4) + y(1)*y(2) &
2092 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2093 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2095 & 0) = (y(2)*y(3)*(y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2096 & + y(2) + y(3) + y(4)))
2099 y(1) = s_cb(i + 3) - s_cb(i + 2)
2100 y(2) = s_cb(i + 2) - s_cb(i + 1)
2101 y(3) = s_cb(i + 1) - s_cb(i)
2102 y(4) = s_cb(i) - s_cb(i - 1)
2104 & 2) = (y(3)*(y(2) + y(3))*(y(1) + y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) &
2105 & + y(4))*(y(1) + y(2) + y(3) + y(4)))
2107 & 1) = (y(3)*y(4)*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2108 & + 6*y(2)*y(3) + 2*y(4)*y(2) + 3*y(3)**2 + 2*y(4)*y(3)))/((y(2) + y(3))*(y(1) &
2109 & + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2111 & 0) = -(y(3)*y(4)*(y(2) + y(3)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2112 & + y(2) + y(3) + y(4)))
2115 y(1) = s_cb(i + 4) - s_cb(i + 3)
2116 y(2) = s_cb(i + 3) - s_cb(i + 2)
2117 y(3) = s_cb(i + 2) - s_cb(i + 1)
2118 y(4) = s_cb(i + 1) - s_cb(i)
2120 & 2) = (y(4)*(y(2)**2 + 4*y(2)*y(3) + 4*y(2)*y(4) + y(1)*y(2) + 3*y(3)**2 &
2121 & + 6*y(3)*y(4) + 2*y(1)*y(3) + 3*y(4)**2 + 2*y(1)*y(4)))/((y(3) + y(4))*(y(2) &
2122 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2124 & 1) = -(y(4)*(y(3) + y(4))*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + 2*y(1)*y(4) &
2125 & + 3*y(2)**2 + 6*y(2)*y(3) + 4*y(2)*y(4) + 3*y(3)**2 + 4*y(3)*y(4) + y(4)**2)) &
2126 & /((y(2) + y(3))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2129 & 0) = (y(4)*(y(3) + y(4))*(y(2) + y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) &
2130 & + y(3))*(y(1) + y(2) + y(3) + y(4)))
2135 y = s_cb(i - 2:i + 1) - s_cb(i - 3:i)
2137 & 0) = (4*y(4)**2*(5*y(1)**2*y(2)**2 + 20*y(1)**2*y(2)*y(3) + 15*y(1)**2*y(2)*y(4) &
2138 & + 20*y(1)**2*y(3)**2 + 30*y(1)**2*y(3)*y(4) + 60*y(1)**2*y(4)**2 + 10*y(1)*y(2) &
2139 & **3 + 60*y(1)*y(2)**2*y(3) + 45*y(1)*y(2)**2*y(4) + 110*y(1)*y(2)*y(3)**2 &
2140 & + 165*y(1)*y(2)*y(3)*y(4) + 260*y(1)*y(2)*y(4)**2 + 60*y(1)*y(3)**3 + 135*y(1) &
2141 & *y(3)**2*y(4) + 400*y(1)*y(3)*y(4)**2 + 225*y(1)*y(4)**3 + 5*y(2)**4 + 40*y(2) &
2142 & **3*y(3) + 30*y(2)**3*y(4) + 110*y(2)**2*y(3)**2 + 165*y(2)**2*y(3)*y(4) &
2143 & + 260*y(2)**2*y(4)**2 + 120*y(2)*y(3)**3 + 270*y(2)*y(3)**2*y(4) + 800*y(2)*y(3) &
2144 & *y(4)**2 + 450*y(2)*y(4)**3 + 45*y(3)**4 + 135*y(3)**3*y(4) + 600*y(3)**2*y(4) &
2145 & **2 + 675*y(3)*y(4)**3 + 996*y(4)**4))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4)) &
2146 & **2*(y(1) + y(2) + y(3) + y(4))**2)
2148 & 1) = -(4*y(4)**2*(10*y(1)**3*y(2)*y(3) + 5*y(1)**3*y(2)*y(4) + 20*y(1)**3*y(3) &
2149 & **2 + 25*y(1)**3*y(3)*y(4) + 105*y(1)**3*y(4)**2 + 40*y(1)**2*y(2)**2*y(3) &
2150 & + 20*y(1)**2*y(2)**2*y(4) + 130*y(1)**2*y(2)*y(3)**2 + 155*y(1)**2*y(2)*y(3)*y(4) &
2151 & + 535*y(1)**2*y(2)*y(4)**2 + 90*y(1)**2*y(3)**3 + 165*y(1)**2*y(3)**2*y(4) &
2152 & + 790*y(1)**2*y(3)*y(4)**2 + 415*y(1)**2*y(4)**3 + 60*y(1)*y(2)**3*y(3) + 30*y(1) &
2153 & *y(2)**3*y(4) + 270*y(1)*y(2)**2*y(3)**2 + 315*y(1)*y(2)**2*y(3)*y(4) + 975*y(1) &
2154 & *y(2)**2*y(4)**2 + 360*y(1)*y(2)*y(3)**3 + 645*y(1)*y(2)*y(3)**2*y(4) + 2850*y(1) &
2155 & *y(2)*y(3)*y(4)**2 + 1460*y(1)*y(2)*y(4)**3 + 150*y(1)*y(3)**4 + 360*y(1)*y(3) &
2156 & **3*y(4) + 2000*y(1)*y(3)**2*y(4)**2 + 2005*y(1)*y(3)*y(4)**3 + 2077*y(1)*y(4) &
2157 & **4 + 30*y(2)**4*y(3) + 15*y(2)**4*y(4) + 180*y(2)**3*y(3)**2 + 210*y(2)**3*y(3) &
2158 & *y(4) + 650*y(2)**3*y(4)**2 + 360*y(2)**2*y(3)**3 + 645*y(2)**2*y(3)**2*y(4) &
2159 & + 2850*y(2)**2*y(3)*y(4)**2 + 1460*y(2)**2*y(4)**3 + 300*y(2)*y(3)**4 + 720*y(2) &
2160 & *y(3)**3*y(4) + 4000*y(2)*y(3)**2*y(4)**2 + 4010*y(2)*y(3)*y(4)**3 + 4154*y(2) &
2161 & *y(4)**4 + 90*y(3)**5 + 270*y(3)**4*y(4) + 1800*y(3)**3*y(4)**2 + 2655*y(3) &
2162 & **2*y(4)**3 + 4464*y(3)*y(4)**4 + 1767*y(4)**5))/(5*(y(2) + y(3))*(y(3) + y(4)) &
2163 & *(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2165 & 2) = (4*y(4)**2*(10*y(2)**3*y(3) + 5*y(2)**3*y(4) + 50*y(2)**2*y(3)**2 + 60*y(2) &
2166 & **2*y(3)*y(4) + 10*y(1)*y(2)**2*y(3) + 215*y(2)**2*y(4)**2 + 5*y(1)*y(2)**2*y(4) &
2167 & + 70*y(2)*y(3)**3 + 130*y(2)*y(3)**2*y(4) + 30*y(1)*y(2)*y(3)**2 + 775*y(2)*y(3) &
2168 & *y(4)**2 + 35*y(1)*y(2)*y(3)*y(4) + 415*y(2)*y(4)**3 + 110*y(1)*y(2)*y(4)**2 &
2169 & + 30*y(3)**4 + 75*y(3)**3*y(4) + 20*y(1)*y(3)**3 + 665*y(3)**2*y(4)**2 + 35*y(1) &
2170 & *y(3)**2*y(4) + 725*y(3)*y(4)**3 + 220*y(1)*y(3)*y(4)**2 + 1767*y(4)**4 &
2171 & + 105*y(1)*y(4)**3))/(5*(y(1) + y(2))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) &
2172 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
2174 & 3) = (4*y(4)**2*(5*y(1)**4*y(3)**2 + 5*y(1)**4*y(3)*y(4) + 50*y(1)**4*y(4)**2 &
2175 & + 30*y(1)**3*y(2)*y(3)**2 + 30*y(1)**3*y(2)*y(3)*y(4) + 300*y(1)**3*y(2)*y(4)**2 &
2176 & + 30*y(1)**3*y(3)**3 + 45*y(1)**3*y(3)**2*y(4) + 415*y(1)**3*y(3)*y(4)**2 &
2177 & + 200*y(1)**3*y(4)**3 + 75*y(1)**2*y(2)**2*y(3)**2 + 75*y(1)**2*y(2)**2*y(3)*y(4) &
2178 & + 750*y(1)**2*y(2)**2*y(4)**2 + 150*y(1)**2*y(2)*y(3)**3 + 225*y(1)**2*y(2)*y(3) &
2179 & **2*y(4) + 2075*y(1)**2*y(2)*y(3)*y(4)**2 + 1000*y(1)**2*y(2)*y(4)**3 + 75*y(1) &
2180 & **2*y(3)**4 + 150*y(1)**2*y(3)**3*y(4) + 1390*y(1)**2*y(3)**2*y(4)**2 + 1315*y(1) &
2181 & **2*y(3)*y(4)**3 + 1081*y(1)**2*y(4)**4 + 90*y(1)*y(2)**3*y(3)**2 + 90*y(1)*y(2) &
2182 & **3*y(3)*y(4) + 900*y(1)*y(2)**3*y(4)**2 + 270*y(1)*y(2)**2*y(3)**3 + 405*y(1) &
2183 & *y(2)**2*y(3)**2*y(4) + 3735*y(1)*y(2)**2*y(3)*y(4)**2 + 1800*y(1)*y(2)**2*y(4) &
2184 & **3 + 270*y(1)*y(2)*y(3)**4 + 540*y(1)*y(2)*y(3)**3*y(4) + 5025*y(1)*y(2)*y(3) &
2185 & **2*y(4)**2 + 4755*y(1)*y(2)*y(3)*y(4)**3 + 4224*y(1)*y(2)*y(4)**4 + 90*y(1)*y(3) &
2186 & **5 + 225*y(1)*y(3)**4*y(4) + 2190*y(1)*y(3)**3*y(4)**2 + 3060*y(1)*y(3)**2*y(4) &
2187 & **3 + 4529*y(1)*y(3)*y(4)**4 + 1762*y(1)*y(4)**5 + 45*y(2)**4*y(3)**2 + 45*y(2) &
2188 & **4*y(3)*y(4) + 450*y(2)**4*y(4)**2 + 180*y(2)**3*y(3)**3 + 270*y(2)**3*y(3) &
2189 & **2*y(4) + 2490*y(2)**3*y(3)*y(4)**2 + 1200*y(2)**3*y(4)**3 + 270*y(2)**2*y(3) &
2190 & **4 + 540*y(2)**2*y(3)**3*y(4) + 5025*y(2)**2*y(3)**2*y(4)**2 + 4755*y(2)**2*y(3) &
2191 & *y(4)**3 + 4224*y(2)**2*y(4)**4 + 180*y(2)*y(3)**5 + 450*y(2)*y(3)**4*y(4) &
2192 & + 4380*y(2)*y(3)**3*y(4)**2 + 6120*y(2)*y(3)**2*y(4)**3 + 9058*y(2)*y(3)*y(4)**4 &
2193 & + 3524*y(2)*y(4)**5 + 45*y(3)**6 + 135*y(3)**5*y(4) + 1395*y(3)**4*y(4)**2 &
2194 & + 2565*y(3)**3*y(4)**3 + 4884*y(3)**2*y(4)**4 + 3624*y(3)*y(4)**5 + 831*y(4)**6)) &
2195 & /(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2196 & + y(3) + y(4))**2)
2198 & 4) = -(4*y(4)**2*(10*y(1)**2*y(2)*y(3)**2 + 10*y(1)**2*y(2)*y(3)*y(4) + 100*y(1) &
2199 & **2*y(2)*y(4)**2 + 10*y(1)**2*y(3)**3 + 15*y(1)**2*y(3)**2*y(4) + 205*y(1) &
2200 & **2*y(3)*y(4)**2 + 100*y(1)**2*y(4)**3 + 30*y(1)*y(2)**2*y(3)**2 + 30*y(1)*y(2) &
2201 & **2*y(3)*y(4) + 300*y(1)*y(2)**2*y(4)**2 + 60*y(1)*y(2)*y(3)**3 + 90*y(1)*y(2) &
2202 & *y(3)**2*y(4) + 1030*y(1)*y(2)*y(3)*y(4)**2 + 500*y(1)*y(2)*y(4)**3 + 30*y(1) &
2203 & *y(3)**4 + 60*y(1)*y(3)**3*y(4) + 835*y(1)*y(3)**2*y(4)**2 + 805*y(1)*y(3)*y(4) &
2204 & **3 + 1762*y(1)*y(4)**4 + 30*y(2)**3*y(3)**2 + 30*y(2)**3*y(3)*y(4) + 300*y(2) &
2205 & **3*y(4)**2 + 90*y(2)**2*y(3)**3 + 135*y(2)**2*y(3)**2*y(4) + 1445*y(2)**2*y(3) &
2206 & *y(4)**2 + 700*y(2)**2*y(4)**3 + 90*y(2)*y(3)**4 + 180*y(2)*y(3)**3*y(4) &
2207 & + 2205*y(2)*y(3)**2*y(4)**2 + 2115*y(2)*y(3)*y(4)**3 + 3624*y(2)*y(4)**4 &
2208 & + 30*y(3)**5 + 75*y(3)**4*y(4) + 1060*y(3)**3*y(4)**2 + 1515*y(3)**2*y(4)**3 &
2209 & + 3824*y(3)*y(4)**4 + 1662*y(4)**5))/(5*(y(1) + y(2))*(y(2) + y(3))*(y(1) + y(2) &
2210 & + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
2212 & 5) = (4*y(4)**2*(5*y(2)**2*y(3)**2 + 5*y(2)**2*y(3)*y(4) + 50*y(2)**2*y(4)**2 &
2213 & + 10*y(2)*y(3)**3 + 15*y(2)*y(3)**2*y(4) + 205*y(2)*y(3)*y(4)**2 + 100*y(2)*y(4) &
2214 & **3 + 5*y(3)**4 + 10*y(3)**3*y(4) + 205*y(3)**2*y(4)**2 + 200*y(3)*y(4)**3 &
2215 & + 831*y(4)**4))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) + y(3) &
2218 y = s_cb(i - 1:i + 2) - s_cb(i - 2:i + 1)
2220 & 0) = (4*y(3)**2*(5*y(1)**2*y(2)**2 + 5*y(1)**2*y(2)*y(3) + 50*y(1)**2*y(3)**2 &
2221 & + 10*y(1)*y(2)**3 + 15*y(1)*y(2)**2*y(3) + 205*y(1)*y(2)*y(3)**2 + 100*y(1)*y(3) &
2222 & **3 + 5*y(2)**4 + 10*y(2)**3*y(3) + 205*y(2)**2*y(3)**2 + 200*y(2)*y(3)**3 &
2223 & + 831*y(3)**4))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) &
2226 & 1) = (4*y(3)**2*(5*y(1)**3*y(2)*y(3) + 10*y(1)**3*y(2)*y(4) - 95*y(1)**3*y(3)**2 &
2227 & + 5*y(1)**3*y(3)*y(4) + 20*y(1)**2*y(2)**2*y(3) + 40*y(1)**2*y(2)**2*y(4) &
2228 & - 465*y(1)**2*y(2)*y(3)**2 + 55*y(1)**2*y(2)*y(3)*y(4) + 10*y(1)**2*y(2)*y(4)**2 &
2229 & - 285*y(1)**2*y(3)**3 + 20*y(1)**2*y(3)**2*y(4) + 5*y(1)**2*y(3)*y(4)**2 &
2230 & + 30*y(1)*y(2)**3*y(3) + 60*y(1)*y(2)**3*y(4) - 825*y(1)*y(2)**2*y(3)**2 &
2231 & + 135*y(1)*y(2)**2*y(3)*y(4) + 30*y(1)*y(2)**2*y(4)**2 - 1040*y(1)*y(2)*y(3)**3 &
2232 & + 100*y(1)*y(2)*y(3)**2*y(4) + 35*y(1)*y(2)*y(3)*y(4)**2 - 1847*y(1)*y(3)**4 &
2233 & + 125*y(1)*y(3)**3*y(4) + 110*y(1)*y(3)**2*y(4)**2 + 15*y(2)**4*y(3) + 30*y(2) &
2234 & **4*y(4) - 550*y(2)**3*y(3)**2 + 90*y(2)**3*y(3)*y(4) + 20*y(2)**3*y(4)**2 &
2235 & - 1040*y(2)**2*y(3)**3 + 100*y(2)**2*y(3)**2*y(4) + 35*y(2)**2*y(3)*y(4)**2 &
2236 & - 3694*y(2)*y(3)**4 + 250*y(2)*y(3)**3*y(4) + 220*y(2)*y(3)**2*y(4)**2 &
2237 & - 3219*y(3)**5 - 1452*y(3)**4*y(4) + 105*y(3)**3*y(4)**2))/(5*(y(2) + y(3))*(y(3) &
2238 & + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4)) &
2241 & 2) = -(4*y(3)**2*(5*y(2)**3*y(3) - 95*y(2)*y(3)**3 - 190*y(2)**2*y(3)**2 &
2242 & + 10*y(2)**3*y(4) + 100*y(3)**3*y(4) - 1562*y(3)**4 - 95*y(1)*y(2)*y(3)**2 &
2243 & + 5*y(1)*y(2)**2*y(3) + 10*y(1)*y(2)**2*y(4) + 100*y(1)*y(3)**2*y(4) + 205*y(2) &
2244 & *y(3)**2*y(4) + 15*y(2)**2*y(3)*y(4) + 10*y(1)*y(2)*y(3)*y(4)))/(5*(y(1) + y(2)) &
2245 & *(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2248 & 3) = (4*y(3)**2*(50*y(1)**4*y(3)**2 + 5*y(1)**4*y(3)*y(4) + 5*y(1)**4*y(4)**2 &
2249 & + 300*y(1)**3*y(2)*y(3)**2 + 30*y(1)**3*y(2)*y(3)*y(4) + 30*y(1)**3*y(2)*y(4)**2 &
2250 & + 200*y(1)**3*y(3)**3 + 25*y(1)**3*y(3)**2*y(4) + 35*y(1)**3*y(3)*y(4)**2 &
2251 & + 10*y(1)**3*y(4)**3 + 750*y(1)**2*y(2)**2*y(3)**2 + 75*y(1)**2*y(2)**2*y(3)*y(4) &
2252 & + 75*y(1)**2*y(2)**2*y(4)**2 + 1000*y(1)**2*y(2)*y(3)**3 + 125*y(1)**2*y(2)*y(3) &
2253 & **2*y(4) + 175*y(1)**2*y(2)*y(3)*y(4)**2 + 50*y(1)**2*y(2)*y(4)**3 + 1081*y(1) &
2254 & **2*y(3)**4 - 50*y(1)**2*y(3)**3*y(4) - 10*y(1)**2*y(3)**2*y(4)**2 + 45*y(1) &
2255 & **2*y(3)*y(4)**3 + 5*y(1)**2*y(4)**4 + 900*y(1)*y(2)**3*y(3)**2 + 90*y(1)*y(2) &
2256 & **3*y(3)*y(4) + 90*y(1)*y(2)**3*y(4)**2 + 1800*y(1)*y(2)**2*y(3)**3 + 225*y(1) &
2257 & *y(2)**2*y(3)**2*y(4) + 315*y(1)*y(2)**2*y(3)*y(4)**2 + 90*y(1)*y(2)**2*y(4)**3 &
2258 & + 4224*y(1)*y(2)*y(3)**4 - 120*y(1)*y(2)*y(3)**3*y(4) + 25*y(1)*y(2)*y(3)**2*y(4) &
2259 & **2 + 165*y(1)*y(2)*y(3)*y(4)**3 + 20*y(1)*y(2)*y(4)**4 + 3324*y(1)*y(3)**5 &
2260 & + 1407*y(1)*y(3)**4*y(4) - 100*y(1)*y(3)**3*y(4)**2 + 70*y(1)*y(3)**2*y(4)**3 &
2261 & + 15*y(1)*y(3)*y(4)**4 + 450*y(2)**4*y(3)**2 + 45*y(2)**4*y(3)*y(4) + 45*y(2) &
2262 & **4*y(4)**2 + 1200*y(2)**3*y(3)**3 + 150*y(2)**3*y(3)**2*y(4) + 210*y(2)**3*y(3) &
2263 & *y(4)**2 + 60*y(2)**3*y(4)**3 + 4224*y(2)**2*y(3)**4 - 120*y(2)**2*y(3)**3*y(4) &
2264 & + 25*y(2)**2*y(3)**2*y(4)**2 + 165*y(2)**2*y(3)*y(4)**3 + 20*y(2)**2*y(4)**4 &
2265 & + 6648*y(2)*y(3)**5 + 2814*y(2)*y(3)**4*y(4) - 200*y(2)*y(3)**3*y(4)**2 &
2266 & + 140*y(2)*y(3)**2*y(4)**3 + 30*y(2)*y(3)*y(4)**4 + 3174*y(3)**6 + 3039*y(3) &
2267 & **5*y(4) + 771*y(3)**4*y(4)**2 + 135*y(3)**3*y(4)**3 + 60*y(3)**2*y(4)**4)) &
2268 & /(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2269 & + y(3) + y(4))**2)
2271 & 4) = -(4*y(3)**2*(100*y(1)**2*y(2)*y(3)**2 + 10*y(1)**2*y(2)*y(3)*y(4) + 10*y(1) &
2272 & **2*y(2)*y(4)**2 - 95*y(1)**2*y(3)**2*y(4) + 5*y(1)**2*y(3)*y(4)**2 + 300*y(1) &
2273 & *y(2)**2*y(3)**2 + 30*y(1)*y(2)**2*y(3)*y(4) + 30*y(1)*y(2)**2*y(4)**2 + 200*y(1) &
2274 & *y(2)*y(3)**3 - 260*y(1)*y(2)*y(3)**2*y(4) + 50*y(1)*y(2)*y(3)*y(4)**2 + 10*y(1) &
2275 & *y(2)*y(4)**3 + 1562*y(1)*y(3)**4 - 190*y(1)*y(3)**3*y(4) + 15*y(1)*y(3)**2*y(4) &
2276 & **2 + 5*y(1)*y(3)*y(4)**3 + 300*y(2)**3*y(3)**2 + 30*y(2)**3*y(3)*y(4) + 30*y(2) &
2277 & **3*y(4)**2 + 400*y(2)**2*y(3)**3 - 235*y(2)**2*y(3)**2*y(4) + 85*y(2)**2*y(3) &
2278 & *y(4)**2 + 20*y(2)**2*y(4)**3 + 3224*y(2)*y(3)**4 - 460*y(2)*y(3)**3*y(4) &
2279 & - 35*y(2)*y(3)**2*y(4)**2 + 25*y(2)*y(3)*y(4)**3 + 3124*y(3)**5 + 1467*y(3) &
2280 & **4*y(4) + 110*y(3)**3*y(4)**2 + 105*y(3)**2*y(4)**3))/(5*(y(1) + y(2))*(y(2) &
2281 & + y(3))*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)) &
2284 & 5) = (4*y(3)**2*(50*y(2)**2*y(3)**2 + 5*y(2)**2*y(3)*y(4) + 5*y(2)**2*y(4)**2 &
2285 & - 95*y(2)*y(3)**2*y(4) + 5*y(2)*y(3)*y(4)**2 + 781*y(3)**4 + 50*y(3)**2*y(4)**2)) &
2286 & /(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) + y(3) + y(4))**2)
2288 y = s_cb(i:i + 3) - s_cb(i - 1:i + 2)
2290 & 0) = (4*y(2)**2*(50*y(1)**2*y(2)**2 + 5*y(1)**2*y(2)*y(3) + 5*y(1)**2*y(3)**2 &
2291 & - 95*y(1)*y(2)**2*y(3) + 5*y(1)*y(2)*y(3)**2 + 781*y(2)**4 + 50*y(2)**2*y(3)**2)) &
2292 & /(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2294 & 1) = -(4*y(2)**2*(105*y(1)**3*y(2)**2 + 25*y(1)**3*y(2)*y(3) + 5*y(1)**3*y(2) &
2295 & *y(4) + 20*y(1)**3*y(3)**2 + 10*y(1)**3*y(3)*y(4) + 110*y(1)**2*y(2)**3 - 35*y(1) &
2296 & **2*y(2)**2*y(3) + 15*y(1)**2*y(2)**2*y(4) + 85*y(1)**2*y(2)*y(3)**2 + 50*y(1) &
2297 & **2*y(2)*y(3)*y(4) + 5*y(1)**2*y(2)*y(4)**2 + 30*y(1)**2*y(3)**3 + 30*y(1) &
2298 & **2*y(3)**2*y(4) + 10*y(1)**2*y(3)*y(4)**2 + 1467*y(1)*y(2)**4 - 460*y(1)*y(2) &
2299 & **3*y(3) - 190*y(1)*y(2)**3*y(4) - 235*y(1)*y(2)**2*y(3)**2 - 260*y(1)*y(2) &
2300 & **2*y(3)*y(4) - 95*y(1)*y(2)**2*y(4)**2 + 30*y(1)*y(2)*y(3)**3 + 30*y(1)*y(2) &
2301 & *y(3)**2*y(4) + 10*y(1)*y(2)*y(3)*y(4)**2 + 3124*y(2)**5 + 3224*y(2)**4*y(3) &
2302 & + 1562*y(2)**4*y(4) + 400*y(2)**3*y(3)**2 + 200*y(2)**3*y(3)*y(4) + 300*y(2) &
2303 & **2*y(3)**3 + 300*y(2)**2*y(3)**2*y(4) + 100*y(2)**2*y(3)*y(4)**2))/(5*(y(2) &
2304 & + y(3))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2305 & + y(3) + y(4))**2)
2307 & 2) = -(4*y(2)**2*(100*y(1)*y(2)**3 - 190*y(2)**2*y(3)**2 + 10*y(1)*y(3)**3 &
2308 & + 5*y(2)*y(3)**3 - 95*y(2)**3*y(3) - 1562*y(2)**4 + 15*y(1)*y(2)*y(3)**2 &
2309 & + 205*y(1)*y(2)**2*y(3) + 100*y(1)*y(2)**2*y(4) + 10*y(1)*y(3)**2*y(4) + 5*y(2) &
2310 & *y(3)**2*y(4) - 95*y(2)**2*y(3)*y(4) + 10*y(1)*y(2)*y(3)*y(4)))/(5*(y(1) + y(2)) &
2311 & *(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2314 & 3) = (4*y(2)**2*(60*y(1)**4*y(2)**2 + 30*y(1)**4*y(2)*y(3) + 15*y(1)**4*y(2)*y(4) &
2315 & + 20*y(1)**4*y(3)**2 + 20*y(1)**4*y(3)*y(4) + 5*y(1)**4*y(4)**2 + 135*y(1) &
2316 & **3*y(2)**3 + 140*y(1)**3*y(2)**2*y(3) + 70*y(1)**3*y(2)**2*y(4) + 165*y(1) &
2317 & **3*y(2)*y(3)**2 + 165*y(1)**3*y(2)*y(3)*y(4) + 45*y(1)**3*y(2)*y(4)**2 + 60*y(1) &
2318 & **3*y(3)**3 + 90*y(1)**3*y(3)**2*y(4) + 50*y(1)**3*y(3)*y(4)**2 + 10*y(1)**3*y(4) &
2319 & **3 + 771*y(1)**2*y(2)**4 - 200*y(1)**2*y(2)**3*y(3) - 100*y(1)**2*y(2)**3*y(4) &
2320 & + 25*y(1)**2*y(2)**2*y(3)**2 + 25*y(1)**2*y(2)**2*y(3)*y(4) - 10*y(1)**2*y(2) &
2321 & **2*y(4)**2 + 210*y(1)**2*y(2)*y(3)**3 + 315*y(1)**2*y(2)*y(3)**2*y(4) + 175*y(1) &
2322 & **2*y(2)*y(3)*y(4)**2 + 35*y(1)**2*y(2)*y(4)**3 + 45*y(1)**2*y(3)**4 + 90*y(1) &
2323 & **2*y(3)**3*y(4) + 75*y(1)**2*y(3)**2*y(4)**2 + 30*y(1)**2*y(3)*y(4)**3 + 5*y(1) &
2324 & **2*y(4)**4 + 3039*y(1)*y(2)**5 + 2814*y(1)*y(2)**4*y(3) + 1407*y(1)*y(2)**4*y(4) &
2325 & - 120*y(1)*y(2)**3*y(3)**2 - 120*y(1)*y(2)**3*y(3)*y(4) - 50*y(1)*y(2)**3*y(4) &
2326 & **2 + 150*y(1)*y(2)**2*y(3)**3 + 225*y(1)*y(2)**2*y(3)**2*y(4) + 125*y(1)*y(2) &
2327 & **2*y(3)*y(4)**2 + 25*y(1)*y(2)**2*y(4)**3 + 45*y(1)*y(2)*y(3)**4 + 90*y(1)*y(2) &
2328 & *y(3)**3*y(4) + 75*y(1)*y(2)*y(3)**2*y(4)**2 + 30*y(1)*y(2)*y(3)*y(4)**3 + 5*y(1) &
2329 & *y(2)*y(4)**4 + 3174*y(2)**6 + 6648*y(2)**5*y(3) + 3324*y(2)**5*y(4) + 4224*y(2) &
2330 & **4*y(3)**2 + 4224*y(2)**4*y(3)*y(4) + 1081*y(2)**4*y(4)**2 + 1200*y(2)**3*y(3) &
2331 & **3 + 1800*y(2)**3*y(3)**2*y(4) + 1000*y(2)**3*y(3)*y(4)**2 + 200*y(2)**3*y(4) &
2332 & **3 + 450*y(2)**2*y(3)**4 + 900*y(2)**2*y(3)**3*y(4) + 750*y(2)**2*y(3)**2*y(4) &
2333 & **2 + 300*y(2)**2*y(3)*y(4)**3 + 50*y(2)**2*y(4)**4))/(5*(y(2) + y(3))**2*(y(1) &
2334 & + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2336 & 4) = (4*y(2)**2*(105*y(1)**2*y(2)**3 + 220*y(1)**2*y(2)**2*y(3) + 110*y(1) &
2337 & **2*y(2)**2*y(4) + 35*y(1)**2*y(2)*y(3)**2 + 35*y(1)**2*y(2)*y(3)*y(4) + 5*y(1) &
2338 & **2*y(2)*y(4)**2 + 20*y(1)**2*y(3)**3 + 30*y(1)**2*y(3)**2*y(4) + 10*y(1)**2*y(3) &
2339 & *y(4)**2 - 1452*y(1)*y(2)**4 + 250*y(1)*y(2)**3*y(3) + 125*y(1)*y(2)**3*y(4) &
2340 & + 100*y(1)*y(2)**2*y(3)**2 + 100*y(1)*y(2)**2*y(3)*y(4) + 20*y(1)*y(2)**2*y(4) &
2341 & **2 + 90*y(1)*y(2)*y(3)**3 + 135*y(1)*y(2)*y(3)**2*y(4) + 55*y(1)*y(2)*y(3)*y(4) &
2342 & **2 + 5*y(1)*y(2)*y(4)**3 + 30*y(1)*y(3)**4 + 60*y(1)*y(3)**3*y(4) + 40*y(1)*y(3) &
2343 & **2*y(4)**2 + 10*y(1)*y(3)*y(4)**3 - 3219*y(2)**5 - 3694*y(2)**4*y(3) - 1847*y(2) &
2344 & **4*y(4) - 1040*y(2)**3*y(3)**2 - 1040*y(2)**3*y(3)*y(4) - 285*y(2)**3*y(4)**2 &
2345 & - 550*y(2)**2*y(3)**3 - 825*y(2)**2*y(3)**2*y(4) - 465*y(2)**2*y(3)*y(4)**2 &
2346 & - 95*y(2)**2*y(4)**3 + 15*y(2)*y(3)**4 + 30*y(2)*y(3)**3*y(4) + 20*y(2)*y(3) &
2347 & **2*y(4)**2 + 5*y(2)*y(3)*y(4)**3))/(5*(y(1) + y(2))*(y(2) + y(3))*(y(1) + y(2) &
2348 & + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
2350 & 5) = (4*y(2)**2*(831*y(2)**4 + 200*y(2)**3*y(3) + 100*y(2)**3*y(4) + 205*y(2) &
2351 & **2*y(3)**2 + 205*y(2)**2*y(3)*y(4) + 50*y(2)**2*y(4)**2 + 10*y(2)*y(3)**3 &
2352 & + 15*y(2)*y(3)**2*y(4) + 5*y(2)*y(3)*y(4)**2 + 5*y(3)**4 + 10*y(3)**3*y(4) &
2353 & + 5*y(3)**2*y(4)**2))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) &
2354 & + y(3) + y(4))**2)
2356 y = s_cb(i + 1:i + 4) - s_cb(i:i + 3)
2358 & 0) = (4*y(1)**2*(831*y(1)**4 + 200*y(1)**3*y(2) + 100*y(1)**3*y(3) + 205*y(1) &
2359 & **2*y(2)**2 + 205*y(1)**2*y(2)*y(3) + 50*y(1)**2*y(3)**2 + 10*y(1)*y(2)**3 &
2360 & + 15*y(1)*y(2)**2*y(3) + 5*y(1)*y(2)*y(3)**2 + 5*y(2)**4 + 10*y(2)**3*y(3) &
2361 & + 5*y(2)**2*y(3)**2))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2362 & + y(3) + y(4))**2)
2364 & 1) = -(4*y(1)**2*(1662*y(1)**5 + 3824*y(1)**4*y(2) + 3624*y(1)**4*y(3) &
2365 & + 1762*y(1)**4*y(4) + 1515*y(1)**3*y(2)**2 + 2115*y(1)**3*y(2)*y(3) + 805*y(1) &
2366 & **3*y(2)*y(4) + 700*y(1)**3*y(3)**2 + 500*y(1)**3*y(3)*y(4) + 100*y(1)**3*y(4) &
2367 & **2 + 1060*y(1)**2*y(2)**3 + 2205*y(1)**2*y(2)**2*y(3) + 835*y(1)**2*y(2)**2*y(4) &
2368 & + 1445*y(1)**2*y(2)*y(3)**2 + 1030*y(1)**2*y(2)*y(3)*y(4) + 205*y(1)**2*y(2)*y(4) &
2369 & **2 + 300*y(1)**2*y(3)**3 + 300*y(1)**2*y(3)**2*y(4) + 100*y(1)**2*y(3)*y(4)**2 &
2370 & + 75*y(1)*y(2)**4 + 180*y(1)*y(2)**3*y(3) + 60*y(1)*y(2)**3*y(4) + 135*y(1)*y(2) &
2371 & **2*y(3)**2 + 90*y(1)*y(2)**2*y(3)*y(4) + 15*y(1)*y(2)**2*y(4)**2 + 30*y(1)*y(2) &
2372 & *y(3)**3 + 30*y(1)*y(2)*y(3)**2*y(4) + 10*y(1)*y(2)*y(3)*y(4)**2 + 30*y(2)**5 &
2373 & + 90*y(2)**4*y(3) + 30*y(2)**4*y(4) + 90*y(2)**3*y(3)**2 + 60*y(2)**3*y(3)*y(4) &
2374 & + 10*y(2)**3*y(4)**2 + 30*y(2)**2*y(3)**3 + 30*y(2)**2*y(3)**2*y(4) + 10*y(2) &
2375 & **2*y(3)*y(4)**2))/(5*(y(2) + y(3))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) &
2376 & + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2378 & 2) = (4*y(1)**2*(1767*y(1)**4 + 725*y(1)**3*y(2) + 415*y(1)**3*y(3) + 105*y(4) &
2379 & *y(1)**3 + 665*y(1)**2*y(2)**2 + 775*y(1)**2*y(2)*y(3) + 220*y(4)*y(1)**2*y(2) &
2380 & + 215*y(1)**2*y(3)**2 + 110*y(4)*y(1)**2*y(3) + 75*y(1)*y(2)**3 + 130*y(1)*y(2) &
2381 & **2*y(3) + 35*y(4)*y(1)*y(2)**2 + 60*y(1)*y(2)*y(3)**2 + 35*y(4)*y(1)*y(2)*y(3) &
2382 & + 5*y(1)*y(3)**3 + 5*y(4)*y(1)*y(3)**2 + 30*y(2)**4 + 70*y(2)**3*y(3) + 20*y(4) &
2383 & *y(2)**3 + 50*y(2)**2*y(3)**2 + 30*y(4)*y(2)**2*y(3) + 10*y(2)*y(3)**3 + 10*y(4) &
2384 & *y(2)*y(3)**2))/(5*(y(1) + y(2))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) &
2385 & + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
2387 & 3) = (4*y(1)**2*(831*y(1)**6 + 3624*y(1)**5*y(2) + 3524*y(1)**5*y(3) + 1762*y(1) &
2388 & **5*y(4) + 4884*y(1)**4*y(2)**2 + 9058*y(1)**4*y(2)*y(3) + 4529*y(1)**4*y(2)*y(4) &
2389 & + 4224*y(1)**4*y(3)**2 + 4224*y(1)**4*y(3)*y(4) + 1081*y(1)**4*y(4)**2 &
2390 & + 2565*y(1)**3*y(2)**3 + 6120*y(1)**3*y(2)**2*y(3) + 3060*y(1)**3*y(2)**2*y(4) &
2391 & + 4755*y(1)**3*y(2)*y(3)**2 + 4755*y(1)**3*y(2)*y(3)*y(4) + 1315*y(1)**3*y(2) &
2392 & *y(4)**2 + 1200*y(1)**3*y(3)**3 + 1800*y(1)**3*y(3)**2*y(4) + 1000*y(1)**3*y(3) &
2393 & *y(4)**2 + 200*y(1)**3*y(4)**3 + 1395*y(1)**2*y(2)**4 + 4380*y(1)**2*y(2)**3*y(3) &
2394 & + 2190*y(1)**2*y(2)**3*y(4) + 5025*y(1)**2*y(2)**2*y(3)**2 + 5025*y(1)**2*y(2) &
2395 & **2*y(3)*y(4) + 1390*y(1)**2*y(2)**2*y(4)**2 + 2490*y(1)**2*y(2)*y(3)**3 &
2396 & + 3735*y(1)**2*y(2)*y(3)**2*y(4) + 2075*y(1)**2*y(2)*y(3)*y(4)**2 + 415*y(1) &
2397 & **2*y(2)*y(4)**3 + 450*y(1)**2*y(3)**4 + 900*y(1)**2*y(3)**3*y(4) + 750*y(1) &
2398 & **2*y(3)**2*y(4)**2 + 300*y(1)**2*y(3)*y(4)**3 + 50*y(1)**2*y(4)**4 + 135*y(1) &
2399 & *y(2)**5 + 450*y(1)*y(2)**4*y(3) + 225*y(1)*y(2)**4*y(4) + 540*y(1)*y(2)**3*y(3) &
2400 & **2 + 540*y(1)*y(2)**3*y(3)*y(4) + 150*y(1)*y(2)**3*y(4)**2 + 270*y(1)*y(2) &
2401 & **2*y(3)**3 + 405*y(1)*y(2)**2*y(3)**2*y(4) + 225*y(1)*y(2)**2*y(3)*y(4)**2 &
2402 & + 45*y(1)*y(2)**2*y(4)**3 + 45*y(1)*y(2)*y(3)**4 + 90*y(1)*y(2)*y(3)**3*y(4) &
2403 & + 75*y(1)*y(2)*y(3)**2*y(4)**2 + 30*y(1)*y(2)*y(3)*y(4)**3 + 5*y(1)*y(2)*y(4)**4 &
2404 & + 45*y(2)**6 + 180*y(2)**5*y(3) + 90*y(2)**5*y(4) + 270*y(2)**4*y(3)**2 &
2405 & + 270*y(2)**4*y(3)*y(4) + 75*y(2)**4*y(4)**2 + 180*y(2)**3*y(3)**3 + 270*y(2) &
2406 & **3*y(3)**2*y(4) + 150*y(2)**3*y(3)*y(4)**2 + 30*y(2)**3*y(4)**3 + 45*y(2) &
2407 & **2*y(3)**4 + 90*y(2)**2*y(3)**3*y(4) + 75*y(2)**2*y(3)**2*y(4)**2 + 30*y(2) &
2408 & **2*y(3)*y(4)**3 + 5*y(2)**2*y(4)**4))/(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3)) &
2409 & **2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2411 & 4) = -(4*y(1)**2*(1767*y(1)**5 + 4464*y(1)**4*y(2) + 4154*y(1)**4*y(3) &
2412 & + 2077*y(1)**4*y(4) + 2655*y(1)**3*y(2)**2 + 4010*y(1)**3*y(2)*y(3) + 2005*y(1) &
2413 & **3*y(2)*y(4) + 1460*y(1)**3*y(3)**2 + 1460*y(1)**3*y(3)*y(4) + 415*y(1)**3*y(4) &
2414 & **2 + 1800*y(1)**2*y(2)**3 + 4000*y(1)**2*y(2)**2*y(3) + 2000*y(1)**2*y(2) &
2415 & **2*y(4) + 2850*y(1)**2*y(2)*y(3)**2 + 2850*y(1)**2*y(2)*y(3)*y(4) + 790*y(1) &
2416 & **2*y(2)*y(4)**2 + 650*y(1)**2*y(3)**3 + 975*y(1)**2*y(3)**2*y(4) + 535*y(1) &
2417 & **2*y(3)*y(4)**2 + 105*y(1)**2*y(4)**3 + 270*y(1)*y(2)**4 + 720*y(1)*y(2)**3*y(3) &
2418 & + 360*y(1)*y(2)**3*y(4) + 645*y(1)*y(2)**2*y(3)**2 + 645*y(1)*y(2)**2*y(3)*y(4) &
2419 & + 165*y(1)*y(2)**2*y(4)**2 + 210*y(1)*y(2)*y(3)**3 + 315*y(1)*y(2)*y(3)**2*y(4) &
2420 & + 155*y(1)*y(2)*y(3)*y(4)**2 + 25*y(1)*y(2)*y(4)**3 + 15*y(1)*y(3)**4 + 30*y(1) &
2421 & *y(3)**3*y(4) + 20*y(1)*y(3)**2*y(4)**2 + 5*y(1)*y(3)*y(4)**3 + 90*y(2)**5 &
2422 & + 300*y(2)**4*y(3) + 150*y(2)**4*y(4) + 360*y(2)**3*y(3)**2 + 360*y(2)**3*y(3) &
2423 & *y(4) + 90*y(2)**3*y(4)**2 + 180*y(2)**2*y(3)**3 + 270*y(2)**2*y(3)**2*y(4) &
2424 & + 130*y(2)**2*y(3)*y(4)**2 + 20*y(2)**2*y(4)**3 + 30*y(2)*y(3)**4 + 60*y(2)*y(3) &
2425 & **3*y(4) + 40*y(2)*y(3)**2*y(4)**2 + 10*y(2)*y(3)*y(4)**3))/(5*(y(1) + y(2)) &
2426 & *(y(2) + y(3))*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2429 & 5) = (4*y(1)**2*(996*y(1)**4 + 675*y(1)**3*y(2) + 450*y(1)**3*y(3) + 225*y(1) &
2430 & **3*y(4) + 600*y(1)**2*y(2)**2 + 800*y(1)**2*y(2)*y(3) + 400*y(1)**2*y(2)*y(4) &
2431 & + 260*y(1)**2*y(3)**2 + 260*y(1)**2*y(3)*y(4) + 60*y(1)**2*y(4)**2 + 135*y(1) &
2432 & *y(2)**3 + 270*y(1)*y(2)**2*y(3) + 135*y(1)*y(2)**2*y(4) + 165*y(1)*y(2)*y(3)**2 &
2433 & + 165*y(1)*y(2)*y(3)*y(4) + 30*y(1)*y(2)*y(4)**2 + 30*y(1)*y(3)**3 + 45*y(1)*y(3) &
2434 & **2*y(4) + 15*y(1)*y(3)*y(4)**2 + 45*y(2)**4 + 120*y(2)**3*y(3) + 60*y(2)**3*y(4) &
2435 & + 110*y(2)**2*y(3)**2 + 110*y(2)**2*y(3)*y(4) + 20*y(2)**2*y(4)**2 + 40*y(2)*y(3) &
2436 & **3 + 60*y(2)*y(3)**2*y(4) + 20*y(2)*y(3)*y(4)**2 + 5*y(3)**4 + 10*y(3)**3*y(4) &
2437 & + 5*y(3)**2*y(4)**2))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) &
2438 & + y(3) + y(4))**2)
2456# 192 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
2458 if (weno_dir == 3)
then
2459 if (weno_order == 3)
then
2460 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
2462 poly_coef_cbr_z(i + 1, 0, 0) = (s_cb(i) - s_cb(i + 1))/(s_cb(i) - s_cb(i + 2))
2463 poly_coef_cbr_z(i + 1, 1, 0) = (s_cb(i) - s_cb(i + 1))/(s_cb(i - 1) - s_cb(i + 1))
2469 d_cbr_z(0, i + 1) = (s_cb(i - 1) - s_cb(i + 1))/(s_cb(i - 1) - s_cb(i + 2))
2470 d_cbl_z(0, i + 1) = (s_cb(i - 1) - s_cb(i))/(s_cb(i - 1) - s_cb(i + 2))
2476 beta_coef_z(i + 1, 0, 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp/(s_cb(i) - s_cb(i + 2))**2._wp
2477 beta_coef_z(i + 1, 1, 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp/(s_cb(i - 1) - s_cb(i + 1))**2._wp
2482 if (null_weights)
then
2483 if (bc_s%beg == bc_riemann_extrap)
then
2488 if (bc_s%end == bc_riemann_extrap)
then
2496 else if (weno_order == 5)
then
2497 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
2500 & 0) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i + 2)))/((s_cb(i) - s_cb(i &
2501 & + 3))*(s_cb(i + 3) - s_cb(i + 1)))
2503 & 0) = ((s_cb(i - 1) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i)))/((s_cb(i - 1) &
2504 & - s_cb(i + 2))*(s_cb(i + 2) - s_cb(i)))
2506 & 1) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i + 2)))/((s_cb(i - 1) &
2507 & - s_cb(i + 1))*(s_cb(i - 1) - s_cb(i + 2)))
2509 & 1) = ((s_cb(i) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))/((s_cb(i - 2) &
2510 & - s_cb(i))*(s_cb(i - 2) - s_cb(i + 1)))
2512 & 0) = ((s_cb(i + 1) - s_cb(i))*(s_cb(i) - s_cb(i + 2)))/((s_cb(i) - s_cb(i + 3)) &
2513 & *(s_cb(i + 3) - s_cb(i + 1)))
2515 & 0) = ((s_cb(i) - s_cb(i - 1))*(s_cb(i) - s_cb(i + 1)))/((s_cb(i - 1) - s_cb(i &
2516 & + 2))*(s_cb(i) - s_cb(i + 2)))
2518 & 1) = ((s_cb(i + 1) - s_cb(i))*(s_cb(i) - s_cb(i + 2)))/((s_cb(i - 1) - s_cb(i &
2519 & + 1))*(s_cb(i - 1) - s_cb(i + 2)))
2521 & 1) = ((s_cb(i - 1) - s_cb(i))*(s_cb(i) - s_cb(i + 1)))/((s_cb(i - 2) - s_cb(i)) &
2522 & *(s_cb(i - 2) - s_cb(i + 1)))
2525 & 1) = ((s_cb(i) - s_cb(i + 2)) + (s_cb(i + 1) - s_cb(i + 3)))/((s_cb(i) - s_cb(i &
2526 & + 2))*(s_cb(i) - s_cb(i + 3)))*((s_cb(i) - s_cb(i + 1)))
2528 & 0) = ((s_cb(i - 2) - s_cb(i + 1)) + (s_cb(i - 1) - s_cb(i + 1)))/((s_cb(i - 1) &
2529 & - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 2)))*((s_cb(i + 1) - s_cb(i)))
2531 & 1) = ((s_cb(i) - s_cb(i + 2)) + (s_cb(i) - s_cb(i + 3)))/((s_cb(i) - s_cb(i + 2)) &
2532 & *(s_cb(i) - s_cb(i + 3)))*((s_cb(i + 1) - s_cb(i)))
2534 & 0) = ((s_cb(i - 2) - s_cb(i)) + (s_cb(i - 1) - s_cb(i + 1)))/((s_cb(i - 2) &
2535 & - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))*((s_cb(i) - s_cb(i + 1)))
2539 & i + 1) = ((s_cb(i - 2) - s_cb(i + 1))*(s_cb(i + 1) - s_cb(i - 1)))/((s_cb(i - 2) &
2540 & - s_cb(i + 3))*(s_cb(i + 3) - s_cb(i - 1)))
2542 & i + 1) = ((s_cb(i + 1) - s_cb(i + 2))*(s_cb(i + 1) - s_cb(i + 3)))/((s_cb(i - 2) &
2543 & - s_cb(i + 2))*(s_cb(i - 2) - s_cb(i + 3)))
2545 & i + 1) = ((s_cb(i - 2) - s_cb(i))*(s_cb(i) - s_cb(i - 1)))/((s_cb(i - 2) - s_cb(i + 3)) &
2546 & *(s_cb(i + 3) - s_cb(i - 1)))
2548 & i + 1) = ((s_cb(i) - s_cb(i + 2))*(s_cb(i) - s_cb(i + 3)))/((s_cb(i - 2) - s_cb(i + 2)) &
2549 & *(s_cb(i - 2) - s_cb(i + 3)))
2556 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2557 & + (s_cb(i + 1) - s_cb(i))*(s_cb(i + 2) - s_cb(i + 1)) + (s_cb(i + 2) - s_cb(i + 1)) &
2558 & **2._wp)/((s_cb(i) - s_cb(i + 3))**2._wp*(s_cb(i + 1) - s_cb(i + 3))**2._wp)
2561 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(19._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2562 & - (s_cb(i + 1) - s_cb(i))*(s_cb(i + 3) - s_cb(i + 1)) + 2._wp*(s_cb(i + 2) - s_cb(i)) &
2563 & *((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1))))/((s_cb(i) - s_cb(i + 2)) &
2564 & *(s_cb(i) - s_cb(i + 3))**2._wp*(s_cb(i + 3) - s_cb(i + 1)))
2567 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2568 & + (s_cb(i + 1) - s_cb(i))*((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1))) &
2569 & + ((s_cb(i + 2) - s_cb(i)) + (s_cb(i + 3) - s_cb(i + 1)))**2._wp)/((s_cb(i) - s_cb(i &
2570 & + 2))**2._wp*(s_cb(i) - s_cb(i + 3))**2._wp)
2573 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2574 & + (s_cb(i) - s_cb(i - 1))**2._wp + (s_cb(i) - s_cb(i - 1))*(s_cb(i + 1) - s_cb(i))) &
2575 & /((s_cb(i - 1) - s_cb(i + 2))**2._wp*(s_cb(i) - s_cb(i + 2))**2._wp)
2578 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*((s_cb(i) - s_cb(i + 1))*((s_cb(i) &
2579 & - s_cb(i - 1)) + 20._wp*(s_cb(i + 1) - s_cb(i))) + (2._wp*(s_cb(i) - s_cb(i - 1)) &
2580 & + (s_cb(i + 1) - s_cb(i)))*(s_cb(i + 2) - s_cb(i)))/((s_cb(i + 1) - s_cb(i - 1)) &
2581 & *(s_cb(i - 1) - s_cb(i + 2))**2._wp*(s_cb(i + 2) - s_cb(i)))
2584 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2585 & + (s_cb(i + 1) - s_cb(i))*(s_cb(i + 2) - s_cb(i + 1)) + (s_cb(i + 2) - s_cb(i + 1)) &
2586 & **2._wp)/((s_cb(i - 1) - s_cb(i + 1))**2._wp*(s_cb(i - 1) - s_cb(i + 2))**2._wp)
2589 & 0) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(12._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2590 & + ((s_cb(i) - s_cb(i - 2)) + (s_cb(i) - s_cb(i - 1)))**2._wp + 3._wp*((s_cb(i) &
2591 & - s_cb(i - 2)) + (s_cb(i) - s_cb(i - 1)))*(s_cb(i + 1) - s_cb(i)))/((s_cb(i - 2) &
2592 & - s_cb(i + 1))**2._wp*(s_cb(i - 1) - s_cb(i + 1))**2._wp)
2595 & 1) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(19._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2596 & + ((s_cb(i) - s_cb(i - 2))*(s_cb(i) - s_cb(i + 1))) + 2._wp*(s_cb(i + 1) - s_cb(i &
2597 & - 1))*((s_cb(i) - s_cb(i - 2)) + (s_cb(i + 1) - s_cb(i - 1))))/((s_cb(i - 2) &
2598 & - s_cb(i))*(s_cb(i - 2) - s_cb(i + 1))**2._wp*(s_cb(i + 1) - s_cb(i - 1)))
2601 & 2) = 4._wp*(s_cb(i) - s_cb(i + 1))**2._wp*(10._wp*(s_cb(i + 1) - s_cb(i))**2._wp &
2602 & + (s_cb(i) - s_cb(i - 1))**2._wp + (s_cb(i) - s_cb(i - 1))*(s_cb(i + 1) - s_cb(i))) &
2603 & /((s_cb(i - 2) - s_cb(i))**2._wp*(s_cb(i - 2) - s_cb(i + 1))**2._wp)
2608 if (null_weights)
then
2609 if (bc_s%beg == bc_riemann_extrap)
then
2616 if (bc_s%end == bc_riemann_extrap)
then
2618 & s - 1)/sum(
d_cbr_z(:,s - 1))
2620 & s - 1)/sum(
d_cbl_z(:,s - 1))
2626 if (.not. teno)
then
2627 do i = is%beg - 1 + weno_polyn, is%end - 1 - weno_polyn
2640 w = s_cb(i - 3:i + 4) - s_cb(i)
2642 & i + 1) = ((w(5) - w(6))*(w(5) - w(7))*(w(5) - w(8)))/((w(1) - w(6))*(w(1) - w(7)) &
2645 & i + 1) = ((w(1) - w(5))*(w(5) - w(7))*(w(5) - w(8))*(w(1)*w(2) - w(1)*w(6) - w(1) &
2646 & *w(7) - w(2)*w(6) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) + w(6)*w(7) + w(6)*w(8) + w(7) &
2647 & *w(8) + w(1)**2 + w(2)**2))/((w(1) - w(6))*(w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7)) &
2650 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(5) - w(8))*(w(1)*w(2) + w(1)*w(3) + w(2) &
2651 & *w(3) - w(1)*w(7) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) - w(3)*w(7) - w(3)*w(8) + w(7) &
2652 & *w(8) + w(7)**2 + w(8)**2))/((w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7))*(w(2) - w(8)) &
2655 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(3) - w(5)))/((w(1) - w(8))*(w(2) - w(8)) &
2662 w(1) = s_cb(i + 4) - s_cb(i)
2663 w(2) = s_cb(i + 3) - s_cb(i)
2664 w(3) = s_cb(i + 2) - s_cb(i)
2665 w(4) = s_cb(i + 1) - s_cb(i)
2666 w(5) = s_cb(i) - s_cb(i)
2667 w(6) = s_cb(i - 1) - s_cb(i)
2668 w(7) = s_cb(i - 2) - s_cb(i)
2669 w(8) = s_cb(i - 3) - s_cb(i)
2671 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(3) - w(5)))/((w(1) - w(8))*(w(2) - w(8)) &
2674 & i + 1) = ((w(1) - w(5))*(w(2) - w(5))*(w(5) - w(8))*(w(1)*w(2) + w(1)*w(3) + w(2) &
2675 & *w(3) - w(1)*w(7) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) - w(3)*w(7) - w(3)*w(8) + w(7) &
2676 & *w(8) + w(7)**2 + w(8)**2))/((w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7))*(w(2) - w(8)) &
2679 & i + 1) = ((w(1) - w(5))*(w(5) - w(7))*(w(5) - w(8))*(w(1)*w(2) - w(1)*w(6) - w(1) &
2680 & *w(7) - w(2)*w(6) - w(1)*w(8) - w(2)*w(7) - w(2)*w(8) + w(6)*w(7) + w(6)*w(8) + w(7) &
2681 & *w(8) + w(1)**2 + w(2)**2))/((w(1) - w(6))*(w(1) - w(7))*(w(1) - w(8))*(w(2) - w(7)) &
2684 & i + 1) = ((w(5) - w(6))*(w(5) - w(7))*(w(5) - w(8)))/((w(1) - w(6))*(w(1) - w(7)) &
2688 y = s_cb(i + 1:i + 4) - s_cb(i:i + 3)
2690 & 0) = (y(1)*y(2)*(y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2691 & + y(2) + y(3) + y(4)))
2693 & 1) = -(y(1)*y(2)*(3*y(2)**2 + 6*y(2)*y(3) + 3*y(2)*y(4) + 2*y(1)*y(2) &
2694 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2695 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2697 & 2) = (y(1)*(y(1)**2 + 3*y(1)*y(2) + 2*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2698 & + 4*y(2)*y(3) + 2*y(4)*y(2) + y(3)**2 + y(4)*y(3)))/((y(1) + y(2))*(y(1) &
2699 & + y(2) + y(3))*(y(1) + y(2) + y(3) + y(4)))
2701 y = s_cb(i:i + 3) - s_cb(i - 1:i + 2)
2703 & 0) = -(y(2)*y(3)*(y(1) + y(2)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2704 & + y(2) + y(3) + y(4)))
2706 & 1) = (y(2)*(y(1) + y(2))*(y(2)**2 + 4*y(2)*y(3) + 2*y(2)*y(4) + y(1)*y(2) &
2707 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2708 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2710 & 2) = (y(2)*y(3)*(y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2711 & + y(2) + y(3) + y(4)))
2713 y = s_cb(i - 1:i + 2) - s_cb(i - 2:i + 1)
2715 & 0) = (y(3)*(y(2) + y(3))*(y(1) + y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) &
2716 & + y(4))*(y(1) + y(2) + y(3) + y(4)))
2718 & 1) = (y(3)*y(4)*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2719 & + 6*y(2)*y(3) + 2*y(4)*y(2) + 3*y(3)**2 + 2*y(4)*y(3)))/((y(2) + y(3))*(y(1) &
2720 & + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2722 & 2) = -(y(3)*y(4)*(y(2) + y(3)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2723 & + y(2) + y(3) + y(4)))
2725 y = s_cb(i - 2:i + 1) - s_cb(i - 3:i)
2727 & 0) = (y(4)*(y(2)**2 + 4*y(2)*y(3) + 4*y(2)*y(4) + y(1)*y(2) + 3*y(3)**2 &
2728 & + 6*y(3)*y(4) + 2*y(1)*y(3) + 3*y(4)**2 + 2*y(1)*y(4)))/((y(3) + y(4))*(y(2) &
2729 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2731 & 1) = -(y(4)*(y(3) + y(4))*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + 2*y(1)*y(4) &
2732 & + 3*y(2)**2 + 6*y(2)*y(3) + 4*y(2)*y(4) + 3*y(3)**2 + 4*y(3)*y(4) + y(4)**2)) &
2733 & /((y(2) + y(3))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2736 & 2) = (y(4)*(y(3) + y(4))*(y(2) + y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) &
2737 & + y(3))*(y(1) + y(2) + y(3) + y(4)))
2740 y(1) = s_cb(i + 1) - s_cb(i)
2741 y(2) = s_cb(i) - s_cb(i - 1)
2742 y(3) = s_cb(i - 1) - s_cb(i - 2)
2743 y(4) = s_cb(i - 2) - s_cb(i - 3)
2745 & 2) = (y(1)*y(2)*(y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2746 & + y(2) + y(3) + y(4)))
2748 & 1) = -(y(1)*y(2)*(3*y(2)**2 + 6*y(2)*y(3) + 3*y(2)*y(4) + 2*y(1)*y(2) &
2749 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2750 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2752 & 0) = (y(1)*(y(1)**2 + 3*y(1)*y(2) + 2*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2753 & + 4*y(2)*y(3) + 2*y(4)*y(2) + y(3)**2 + y(4)*y(3)))/((y(1) + y(2))*(y(1) &
2754 & + y(2) + y(3))*(y(1) + y(2) + y(3) + y(4)))
2757 y(1) = s_cb(i + 2) - s_cb(i + 1)
2758 y(2) = s_cb(i + 1) - s_cb(i)
2759 y(3) = s_cb(i) - s_cb(i - 1)
2760 y(4) = s_cb(i - 1) - s_cb(i - 2)
2762 & 2) = -(y(2)*y(3)*(y(1) + y(2)))/((y(3) + y(4))*(y(2) + y(3) + y(4))*(y(1) &
2763 & + y(2) + y(3) + y(4)))
2765 & 1) = (y(2)*(y(1) + y(2))*(y(2)**2 + 4*y(2)*y(3) + 2*y(2)*y(4) + y(1)*y(2) &
2766 & + 3*y(3)**2 + 3*y(3)*y(4) + 2*y(1)*y(3) + y(4)**2 + y(1)*y(4)))/((y(2) + y(3) &
2767 & )*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2769 & 0) = (y(2)*y(3)*(y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2770 & + y(2) + y(3) + y(4)))
2773 y(1) = s_cb(i + 3) - s_cb(i + 2)
2774 y(2) = s_cb(i + 2) - s_cb(i + 1)
2775 y(3) = s_cb(i + 1) - s_cb(i)
2776 y(4) = s_cb(i) - s_cb(i - 1)
2778 & 2) = (y(3)*(y(2) + y(3))*(y(1) + y(2) + y(3)))/((y(3) + y(4))*(y(2) + y(3) &
2779 & + y(4))*(y(1) + y(2) + y(3) + y(4)))
2781 & 1) = (y(3)*y(4)*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + y(4)*y(1) + 3*y(2)**2 &
2782 & + 6*y(2)*y(3) + 2*y(4)*y(2) + 3*y(3)**2 + 2*y(4)*y(3)))/((y(2) + y(3))*(y(1) &
2783 & + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2785 & 0) = -(y(3)*y(4)*(y(2) + y(3)))/((y(1) + y(2))*(y(1) + y(2) + y(3))*(y(1) &
2786 & + y(2) + y(3) + y(4)))
2789 y(1) = s_cb(i + 4) - s_cb(i + 3)
2790 y(2) = s_cb(i + 3) - s_cb(i + 2)
2791 y(3) = s_cb(i + 2) - s_cb(i + 1)
2792 y(4) = s_cb(i + 1) - s_cb(i)
2794 & 2) = (y(4)*(y(2)**2 + 4*y(2)*y(3) + 4*y(2)*y(4) + y(1)*y(2) + 3*y(3)**2 &
2795 & + 6*y(3)*y(4) + 2*y(1)*y(3) + 3*y(4)**2 + 2*y(1)*y(4)))/((y(3) + y(4))*(y(2) &
2796 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)))
2798 & 1) = -(y(4)*(y(3) + y(4))*(y(1)**2 + 3*y(1)*y(2) + 3*y(1)*y(3) + 2*y(1)*y(4) &
2799 & + 3*y(2)**2 + 6*y(2)*y(3) + 4*y(2)*y(4) + 3*y(3)**2 + 4*y(3)*y(4) + y(4)**2)) &
2800 & /((y(2) + y(3))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2803 & 0) = (y(4)*(y(3) + y(4))*(y(2) + y(3) + y(4)))/((y(1) + y(2))*(y(1) + y(2) &
2804 & + y(3))*(y(1) + y(2) + y(3) + y(4)))
2809 y = s_cb(i - 2:i + 1) - s_cb(i - 3:i)
2811 & 0) = (4*y(4)**2*(5*y(1)**2*y(2)**2 + 20*y(1)**2*y(2)*y(3) + 15*y(1)**2*y(2)*y(4) &
2812 & + 20*y(1)**2*y(3)**2 + 30*y(1)**2*y(3)*y(4) + 60*y(1)**2*y(4)**2 + 10*y(1)*y(2) &
2813 & **3 + 60*y(1)*y(2)**2*y(3) + 45*y(1)*y(2)**2*y(4) + 110*y(1)*y(2)*y(3)**2 &
2814 & + 165*y(1)*y(2)*y(3)*y(4) + 260*y(1)*y(2)*y(4)**2 + 60*y(1)*y(3)**3 + 135*y(1) &
2815 & *y(3)**2*y(4) + 400*y(1)*y(3)*y(4)**2 + 225*y(1)*y(4)**3 + 5*y(2)**4 + 40*y(2) &
2816 & **3*y(3) + 30*y(2)**3*y(4) + 110*y(2)**2*y(3)**2 + 165*y(2)**2*y(3)*y(4) &
2817 & + 260*y(2)**2*y(4)**2 + 120*y(2)*y(3)**3 + 270*y(2)*y(3)**2*y(4) + 800*y(2)*y(3) &
2818 & *y(4)**2 + 450*y(2)*y(4)**3 + 45*y(3)**4 + 135*y(3)**3*y(4) + 600*y(3)**2*y(4) &
2819 & **2 + 675*y(3)*y(4)**3 + 996*y(4)**4))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4)) &
2820 & **2*(y(1) + y(2) + y(3) + y(4))**2)
2822 & 1) = -(4*y(4)**2*(10*y(1)**3*y(2)*y(3) + 5*y(1)**3*y(2)*y(4) + 20*y(1)**3*y(3) &
2823 & **2 + 25*y(1)**3*y(3)*y(4) + 105*y(1)**3*y(4)**2 + 40*y(1)**2*y(2)**2*y(3) &
2824 & + 20*y(1)**2*y(2)**2*y(4) + 130*y(1)**2*y(2)*y(3)**2 + 155*y(1)**2*y(2)*y(3)*y(4) &
2825 & + 535*y(1)**2*y(2)*y(4)**2 + 90*y(1)**2*y(3)**3 + 165*y(1)**2*y(3)**2*y(4) &
2826 & + 790*y(1)**2*y(3)*y(4)**2 + 415*y(1)**2*y(4)**3 + 60*y(1)*y(2)**3*y(3) + 30*y(1) &
2827 & *y(2)**3*y(4) + 270*y(1)*y(2)**2*y(3)**2 + 315*y(1)*y(2)**2*y(3)*y(4) + 975*y(1) &
2828 & *y(2)**2*y(4)**2 + 360*y(1)*y(2)*y(3)**3 + 645*y(1)*y(2)*y(3)**2*y(4) + 2850*y(1) &
2829 & *y(2)*y(3)*y(4)**2 + 1460*y(1)*y(2)*y(4)**3 + 150*y(1)*y(3)**4 + 360*y(1)*y(3) &
2830 & **3*y(4) + 2000*y(1)*y(3)**2*y(4)**2 + 2005*y(1)*y(3)*y(4)**3 + 2077*y(1)*y(4) &
2831 & **4 + 30*y(2)**4*y(3) + 15*y(2)**4*y(4) + 180*y(2)**3*y(3)**2 + 210*y(2)**3*y(3) &
2832 & *y(4) + 650*y(2)**3*y(4)**2 + 360*y(2)**2*y(3)**3 + 645*y(2)**2*y(3)**2*y(4) &
2833 & + 2850*y(2)**2*y(3)*y(4)**2 + 1460*y(2)**2*y(4)**3 + 300*y(2)*y(3)**4 + 720*y(2) &
2834 & *y(3)**3*y(4) + 4000*y(2)*y(3)**2*y(4)**2 + 4010*y(2)*y(3)*y(4)**3 + 4154*y(2) &
2835 & *y(4)**4 + 90*y(3)**5 + 270*y(3)**4*y(4) + 1800*y(3)**3*y(4)**2 + 2655*y(3) &
2836 & **2*y(4)**3 + 4464*y(3)*y(4)**4 + 1767*y(4)**5))/(5*(y(2) + y(3))*(y(3) + y(4)) &
2837 & *(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2839 & 2) = (4*y(4)**2*(10*y(2)**3*y(3) + 5*y(2)**3*y(4) + 50*y(2)**2*y(3)**2 + 60*y(2) &
2840 & **2*y(3)*y(4) + 10*y(1)*y(2)**2*y(3) + 215*y(2)**2*y(4)**2 + 5*y(1)*y(2)**2*y(4) &
2841 & + 70*y(2)*y(3)**3 + 130*y(2)*y(3)**2*y(4) + 30*y(1)*y(2)*y(3)**2 + 775*y(2)*y(3) &
2842 & *y(4)**2 + 35*y(1)*y(2)*y(3)*y(4) + 415*y(2)*y(4)**3 + 110*y(1)*y(2)*y(4)**2 &
2843 & + 30*y(3)**4 + 75*y(3)**3*y(4) + 20*y(1)*y(3)**3 + 665*y(3)**2*y(4)**2 + 35*y(1) &
2844 & *y(3)**2*y(4) + 725*y(3)*y(4)**3 + 220*y(1)*y(3)*y(4)**2 + 1767*y(4)**4 &
2845 & + 105*y(1)*y(4)**3))/(5*(y(1) + y(2))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) &
2846 & + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
2848 & 3) = (4*y(4)**2*(5*y(1)**4*y(3)**2 + 5*y(1)**4*y(3)*y(4) + 50*y(1)**4*y(4)**2 &
2849 & + 30*y(1)**3*y(2)*y(3)**2 + 30*y(1)**3*y(2)*y(3)*y(4) + 300*y(1)**3*y(2)*y(4)**2 &
2850 & + 30*y(1)**3*y(3)**3 + 45*y(1)**3*y(3)**2*y(4) + 415*y(1)**3*y(3)*y(4)**2 &
2851 & + 200*y(1)**3*y(4)**3 + 75*y(1)**2*y(2)**2*y(3)**2 + 75*y(1)**2*y(2)**2*y(3)*y(4) &
2852 & + 750*y(1)**2*y(2)**2*y(4)**2 + 150*y(1)**2*y(2)*y(3)**3 + 225*y(1)**2*y(2)*y(3) &
2853 & **2*y(4) + 2075*y(1)**2*y(2)*y(3)*y(4)**2 + 1000*y(1)**2*y(2)*y(4)**3 + 75*y(1) &
2854 & **2*y(3)**4 + 150*y(1)**2*y(3)**3*y(4) + 1390*y(1)**2*y(3)**2*y(4)**2 + 1315*y(1) &
2855 & **2*y(3)*y(4)**3 + 1081*y(1)**2*y(4)**4 + 90*y(1)*y(2)**3*y(3)**2 + 90*y(1)*y(2) &
2856 & **3*y(3)*y(4) + 900*y(1)*y(2)**3*y(4)**2 + 270*y(1)*y(2)**2*y(3)**3 + 405*y(1) &
2857 & *y(2)**2*y(3)**2*y(4) + 3735*y(1)*y(2)**2*y(3)*y(4)**2 + 1800*y(1)*y(2)**2*y(4) &
2858 & **3 + 270*y(1)*y(2)*y(3)**4 + 540*y(1)*y(2)*y(3)**3*y(4) + 5025*y(1)*y(2)*y(3) &
2859 & **2*y(4)**2 + 4755*y(1)*y(2)*y(3)*y(4)**3 + 4224*y(1)*y(2)*y(4)**4 + 90*y(1)*y(3) &
2860 & **5 + 225*y(1)*y(3)**4*y(4) + 2190*y(1)*y(3)**3*y(4)**2 + 3060*y(1)*y(3)**2*y(4) &
2861 & **3 + 4529*y(1)*y(3)*y(4)**4 + 1762*y(1)*y(4)**5 + 45*y(2)**4*y(3)**2 + 45*y(2) &
2862 & **4*y(3)*y(4) + 450*y(2)**4*y(4)**2 + 180*y(2)**3*y(3)**3 + 270*y(2)**3*y(3) &
2863 & **2*y(4) + 2490*y(2)**3*y(3)*y(4)**2 + 1200*y(2)**3*y(4)**3 + 270*y(2)**2*y(3) &
2864 & **4 + 540*y(2)**2*y(3)**3*y(4) + 5025*y(2)**2*y(3)**2*y(4)**2 + 4755*y(2)**2*y(3) &
2865 & *y(4)**3 + 4224*y(2)**2*y(4)**4 + 180*y(2)*y(3)**5 + 450*y(2)*y(3)**4*y(4) &
2866 & + 4380*y(2)*y(3)**3*y(4)**2 + 6120*y(2)*y(3)**2*y(4)**3 + 9058*y(2)*y(3)*y(4)**4 &
2867 & + 3524*y(2)*y(4)**5 + 45*y(3)**6 + 135*y(3)**5*y(4) + 1395*y(3)**4*y(4)**2 &
2868 & + 2565*y(3)**3*y(4)**3 + 4884*y(3)**2*y(4)**4 + 3624*y(3)*y(4)**5 + 831*y(4)**6)) &
2869 & /(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2870 & + y(3) + y(4))**2)
2872 & 4) = -(4*y(4)**2*(10*y(1)**2*y(2)*y(3)**2 + 10*y(1)**2*y(2)*y(3)*y(4) + 100*y(1) &
2873 & **2*y(2)*y(4)**2 + 10*y(1)**2*y(3)**3 + 15*y(1)**2*y(3)**2*y(4) + 205*y(1) &
2874 & **2*y(3)*y(4)**2 + 100*y(1)**2*y(4)**3 + 30*y(1)*y(2)**2*y(3)**2 + 30*y(1)*y(2) &
2875 & **2*y(3)*y(4) + 300*y(1)*y(2)**2*y(4)**2 + 60*y(1)*y(2)*y(3)**3 + 90*y(1)*y(2) &
2876 & *y(3)**2*y(4) + 1030*y(1)*y(2)*y(3)*y(4)**2 + 500*y(1)*y(2)*y(4)**3 + 30*y(1) &
2877 & *y(3)**4 + 60*y(1)*y(3)**3*y(4) + 835*y(1)*y(3)**2*y(4)**2 + 805*y(1)*y(3)*y(4) &
2878 & **3 + 1762*y(1)*y(4)**4 + 30*y(2)**3*y(3)**2 + 30*y(2)**3*y(3)*y(4) + 300*y(2) &
2879 & **3*y(4)**2 + 90*y(2)**2*y(3)**3 + 135*y(2)**2*y(3)**2*y(4) + 1445*y(2)**2*y(3) &
2880 & *y(4)**2 + 700*y(2)**2*y(4)**3 + 90*y(2)*y(3)**4 + 180*y(2)*y(3)**3*y(4) &
2881 & + 2205*y(2)*y(3)**2*y(4)**2 + 2115*y(2)*y(3)*y(4)**3 + 3624*y(2)*y(4)**4 &
2882 & + 30*y(3)**5 + 75*y(3)**4*y(4) + 1060*y(3)**3*y(4)**2 + 1515*y(3)**2*y(4)**3 &
2883 & + 3824*y(3)*y(4)**4 + 1662*y(4)**5))/(5*(y(1) + y(2))*(y(2) + y(3))*(y(1) + y(2) &
2884 & + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
2886 & 5) = (4*y(4)**2*(5*y(2)**2*y(3)**2 + 5*y(2)**2*y(3)*y(4) + 50*y(2)**2*y(4)**2 &
2887 & + 10*y(2)*y(3)**3 + 15*y(2)*y(3)**2*y(4) + 205*y(2)*y(3)*y(4)**2 + 100*y(2)*y(4) &
2888 & **3 + 5*y(3)**4 + 10*y(3)**3*y(4) + 205*y(3)**2*y(4)**2 + 200*y(3)*y(4)**3 &
2889 & + 831*y(4)**4))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) + y(3) &
2892 y = s_cb(i - 1:i + 2) - s_cb(i - 2:i + 1)
2894 & 0) = (4*y(3)**2*(5*y(1)**2*y(2)**2 + 5*y(1)**2*y(2)*y(3) + 50*y(1)**2*y(3)**2 &
2895 & + 10*y(1)*y(2)**3 + 15*y(1)*y(2)**2*y(3) + 205*y(1)*y(2)*y(3)**2 + 100*y(1)*y(3) &
2896 & **3 + 5*y(2)**4 + 10*y(2)**3*y(3) + 205*y(2)**2*y(3)**2 + 200*y(2)*y(3)**3 &
2897 & + 831*y(3)**4))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) &
2900 & 1) = (4*y(3)**2*(5*y(1)**3*y(2)*y(3) + 10*y(1)**3*y(2)*y(4) - 95*y(1)**3*y(3)**2 &
2901 & + 5*y(1)**3*y(3)*y(4) + 20*y(1)**2*y(2)**2*y(3) + 40*y(1)**2*y(2)**2*y(4) &
2902 & - 465*y(1)**2*y(2)*y(3)**2 + 55*y(1)**2*y(2)*y(3)*y(4) + 10*y(1)**2*y(2)*y(4)**2 &
2903 & - 285*y(1)**2*y(3)**3 + 20*y(1)**2*y(3)**2*y(4) + 5*y(1)**2*y(3)*y(4)**2 &
2904 & + 30*y(1)*y(2)**3*y(3) + 60*y(1)*y(2)**3*y(4) - 825*y(1)*y(2)**2*y(3)**2 &
2905 & + 135*y(1)*y(2)**2*y(3)*y(4) + 30*y(1)*y(2)**2*y(4)**2 - 1040*y(1)*y(2)*y(3)**3 &
2906 & + 100*y(1)*y(2)*y(3)**2*y(4) + 35*y(1)*y(2)*y(3)*y(4)**2 - 1847*y(1)*y(3)**4 &
2907 & + 125*y(1)*y(3)**3*y(4) + 110*y(1)*y(3)**2*y(4)**2 + 15*y(2)**4*y(3) + 30*y(2) &
2908 & **4*y(4) - 550*y(2)**3*y(3)**2 + 90*y(2)**3*y(3)*y(4) + 20*y(2)**3*y(4)**2 &
2909 & - 1040*y(2)**2*y(3)**3 + 100*y(2)**2*y(3)**2*y(4) + 35*y(2)**2*y(3)*y(4)**2 &
2910 & - 3694*y(2)*y(3)**4 + 250*y(2)*y(3)**3*y(4) + 220*y(2)*y(3)**2*y(4)**2 &
2911 & - 3219*y(3)**5 - 1452*y(3)**4*y(4) + 105*y(3)**3*y(4)**2))/(5*(y(2) + y(3))*(y(3) &
2912 & + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4)) &
2915 & 2) = -(4*y(3)**2*(5*y(2)**3*y(3) - 95*y(2)*y(3)**3 - 190*y(2)**2*y(3)**2 &
2916 & + 10*y(2)**3*y(4) + 100*y(3)**3*y(4) - 1562*y(3)**4 - 95*y(1)*y(2)*y(3)**2 &
2917 & + 5*y(1)*y(2)**2*y(3) + 10*y(1)*y(2)**2*y(4) + 100*y(1)*y(3)**2*y(4) + 205*y(2) &
2918 & *y(3)**2*y(4) + 15*y(2)**2*y(3)*y(4) + 10*y(1)*y(2)*y(3)*y(4)))/(5*(y(1) + y(2)) &
2919 & *(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2922 & 3) = (4*y(3)**2*(50*y(1)**4*y(3)**2 + 5*y(1)**4*y(3)*y(4) + 5*y(1)**4*y(4)**2 &
2923 & + 300*y(1)**3*y(2)*y(3)**2 + 30*y(1)**3*y(2)*y(3)*y(4) + 30*y(1)**3*y(2)*y(4)**2 &
2924 & + 200*y(1)**3*y(3)**3 + 25*y(1)**3*y(3)**2*y(4) + 35*y(1)**3*y(3)*y(4)**2 &
2925 & + 10*y(1)**3*y(4)**3 + 750*y(1)**2*y(2)**2*y(3)**2 + 75*y(1)**2*y(2)**2*y(3)*y(4) &
2926 & + 75*y(1)**2*y(2)**2*y(4)**2 + 1000*y(1)**2*y(2)*y(3)**3 + 125*y(1)**2*y(2)*y(3) &
2927 & **2*y(4) + 175*y(1)**2*y(2)*y(3)*y(4)**2 + 50*y(1)**2*y(2)*y(4)**3 + 1081*y(1) &
2928 & **2*y(3)**4 - 50*y(1)**2*y(3)**3*y(4) - 10*y(1)**2*y(3)**2*y(4)**2 + 45*y(1) &
2929 & **2*y(3)*y(4)**3 + 5*y(1)**2*y(4)**4 + 900*y(1)*y(2)**3*y(3)**2 + 90*y(1)*y(2) &
2930 & **3*y(3)*y(4) + 90*y(1)*y(2)**3*y(4)**2 + 1800*y(1)*y(2)**2*y(3)**3 + 225*y(1) &
2931 & *y(2)**2*y(3)**2*y(4) + 315*y(1)*y(2)**2*y(3)*y(4)**2 + 90*y(1)*y(2)**2*y(4)**3 &
2932 & + 4224*y(1)*y(2)*y(3)**4 - 120*y(1)*y(2)*y(3)**3*y(4) + 25*y(1)*y(2)*y(3)**2*y(4) &
2933 & **2 + 165*y(1)*y(2)*y(3)*y(4)**3 + 20*y(1)*y(2)*y(4)**4 + 3324*y(1)*y(3)**5 &
2934 & + 1407*y(1)*y(3)**4*y(4) - 100*y(1)*y(3)**3*y(4)**2 + 70*y(1)*y(3)**2*y(4)**3 &
2935 & + 15*y(1)*y(3)*y(4)**4 + 450*y(2)**4*y(3)**2 + 45*y(2)**4*y(3)*y(4) + 45*y(2) &
2936 & **4*y(4)**2 + 1200*y(2)**3*y(3)**3 + 150*y(2)**3*y(3)**2*y(4) + 210*y(2)**3*y(3) &
2937 & *y(4)**2 + 60*y(2)**3*y(4)**3 + 4224*y(2)**2*y(3)**4 - 120*y(2)**2*y(3)**3*y(4) &
2938 & + 25*y(2)**2*y(3)**2*y(4)**2 + 165*y(2)**2*y(3)*y(4)**3 + 20*y(2)**2*y(4)**4 &
2939 & + 6648*y(2)*y(3)**5 + 2814*y(2)*y(3)**4*y(4) - 200*y(2)*y(3)**3*y(4)**2 &
2940 & + 140*y(2)*y(3)**2*y(4)**3 + 30*y(2)*y(3)*y(4)**4 + 3174*y(3)**6 + 3039*y(3) &
2941 & **5*y(4) + 771*y(3)**4*y(4)**2 + 135*y(3)**3*y(4)**3 + 60*y(3)**2*y(4)**4)) &
2942 & /(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2943 & + y(3) + y(4))**2)
2945 & 4) = -(4*y(3)**2*(100*y(1)**2*y(2)*y(3)**2 + 10*y(1)**2*y(2)*y(3)*y(4) + 10*y(1) &
2946 & **2*y(2)*y(4)**2 - 95*y(1)**2*y(3)**2*y(4) + 5*y(1)**2*y(3)*y(4)**2 + 300*y(1) &
2947 & *y(2)**2*y(3)**2 + 30*y(1)*y(2)**2*y(3)*y(4) + 30*y(1)*y(2)**2*y(4)**2 + 200*y(1) &
2948 & *y(2)*y(3)**3 - 260*y(1)*y(2)*y(3)**2*y(4) + 50*y(1)*y(2)*y(3)*y(4)**2 + 10*y(1) &
2949 & *y(2)*y(4)**3 + 1562*y(1)*y(3)**4 - 190*y(1)*y(3)**3*y(4) + 15*y(1)*y(3)**2*y(4) &
2950 & **2 + 5*y(1)*y(3)*y(4)**3 + 300*y(2)**3*y(3)**2 + 30*y(2)**3*y(3)*y(4) + 30*y(2) &
2951 & **3*y(4)**2 + 400*y(2)**2*y(3)**3 - 235*y(2)**2*y(3)**2*y(4) + 85*y(2)**2*y(3) &
2952 & *y(4)**2 + 20*y(2)**2*y(4)**3 + 3224*y(2)*y(3)**4 - 460*y(2)*y(3)**3*y(4) &
2953 & - 35*y(2)*y(3)**2*y(4)**2 + 25*y(2)*y(3)*y(4)**3 + 3124*y(3)**5 + 1467*y(3) &
2954 & **4*y(4) + 110*y(3)**3*y(4)**2 + 105*y(3)**2*y(4)**3))/(5*(y(1) + y(2))*(y(2) &
2955 & + y(3))*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4)) &
2958 & 5) = (4*y(3)**2*(50*y(2)**2*y(3)**2 + 5*y(2)**2*y(3)*y(4) + 5*y(2)**2*y(4)**2 &
2959 & - 95*y(2)*y(3)**2*y(4) + 5*y(2)*y(3)*y(4)**2 + 781*y(3)**4 + 50*y(3)**2*y(4)**2)) &
2960 & /(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) + y(3) + y(4))**2)
2962 y = s_cb(i:i + 3) - s_cb(i - 1:i + 2)
2964 & 0) = (4*y(2)**2*(50*y(1)**2*y(2)**2 + 5*y(1)**2*y(2)*y(3) + 5*y(1)**2*y(3)**2 &
2965 & - 95*y(1)*y(2)**2*y(3) + 5*y(1)*y(2)*y(3)**2 + 781*y(2)**4 + 50*y(2)**2*y(3)**2)) &
2966 & /(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
2968 & 1) = -(4*y(2)**2*(105*y(1)**3*y(2)**2 + 25*y(1)**3*y(2)*y(3) + 5*y(1)**3*y(2) &
2969 & *y(4) + 20*y(1)**3*y(3)**2 + 10*y(1)**3*y(3)*y(4) + 110*y(1)**2*y(2)**3 - 35*y(1) &
2970 & **2*y(2)**2*y(3) + 15*y(1)**2*y(2)**2*y(4) + 85*y(1)**2*y(2)*y(3)**2 + 50*y(1) &
2971 & **2*y(2)*y(3)*y(4) + 5*y(1)**2*y(2)*y(4)**2 + 30*y(1)**2*y(3)**3 + 30*y(1) &
2972 & **2*y(3)**2*y(4) + 10*y(1)**2*y(3)*y(4)**2 + 1467*y(1)*y(2)**4 - 460*y(1)*y(2) &
2973 & **3*y(3) - 190*y(1)*y(2)**3*y(4) - 235*y(1)*y(2)**2*y(3)**2 - 260*y(1)*y(2) &
2974 & **2*y(3)*y(4) - 95*y(1)*y(2)**2*y(4)**2 + 30*y(1)*y(2)*y(3)**3 + 30*y(1)*y(2) &
2975 & *y(3)**2*y(4) + 10*y(1)*y(2)*y(3)*y(4)**2 + 3124*y(2)**5 + 3224*y(2)**4*y(3) &
2976 & + 1562*y(2)**4*y(4) + 400*y(2)**3*y(3)**2 + 200*y(2)**3*y(3)*y(4) + 300*y(2) &
2977 & **2*y(3)**3 + 300*y(2)**2*y(3)**2*y(4) + 100*y(2)**2*y(3)*y(4)**2))/(5*(y(2) &
2978 & + y(3))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
2979 & + y(3) + y(4))**2)
2981 & 2) = -(4*y(2)**2*(100*y(1)*y(2)**3 - 190*y(2)**2*y(3)**2 + 10*y(1)*y(3)**3 &
2982 & + 5*y(2)*y(3)**3 - 95*y(2)**3*y(3) - 1562*y(2)**4 + 15*y(1)*y(2)*y(3)**2 &
2983 & + 205*y(1)*y(2)**2*y(3) + 100*y(1)*y(2)**2*y(4) + 10*y(1)*y(3)**2*y(4) + 5*y(2) &
2984 & *y(3)**2*y(4) - 95*y(2)**2*y(3)*y(4) + 10*y(1)*y(2)*y(3)*y(4)))/(5*(y(1) + y(2)) &
2985 & *(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
2988 & 3) = (4*y(2)**2*(60*y(1)**4*y(2)**2 + 30*y(1)**4*y(2)*y(3) + 15*y(1)**4*y(2)*y(4) &
2989 & + 20*y(1)**4*y(3)**2 + 20*y(1)**4*y(3)*y(4) + 5*y(1)**4*y(4)**2 + 135*y(1) &
2990 & **3*y(2)**3 + 140*y(1)**3*y(2)**2*y(3) + 70*y(1)**3*y(2)**2*y(4) + 165*y(1) &
2991 & **3*y(2)*y(3)**2 + 165*y(1)**3*y(2)*y(3)*y(4) + 45*y(1)**3*y(2)*y(4)**2 + 60*y(1) &
2992 & **3*y(3)**3 + 90*y(1)**3*y(3)**2*y(4) + 50*y(1)**3*y(3)*y(4)**2 + 10*y(1)**3*y(4) &
2993 & **3 + 771*y(1)**2*y(2)**4 - 200*y(1)**2*y(2)**3*y(3) - 100*y(1)**2*y(2)**3*y(4) &
2994 & + 25*y(1)**2*y(2)**2*y(3)**2 + 25*y(1)**2*y(2)**2*y(3)*y(4) - 10*y(1)**2*y(2) &
2995 & **2*y(4)**2 + 210*y(1)**2*y(2)*y(3)**3 + 315*y(1)**2*y(2)*y(3)**2*y(4) + 175*y(1) &
2996 & **2*y(2)*y(3)*y(4)**2 + 35*y(1)**2*y(2)*y(4)**3 + 45*y(1)**2*y(3)**4 + 90*y(1) &
2997 & **2*y(3)**3*y(4) + 75*y(1)**2*y(3)**2*y(4)**2 + 30*y(1)**2*y(3)*y(4)**3 + 5*y(1) &
2998 & **2*y(4)**4 + 3039*y(1)*y(2)**5 + 2814*y(1)*y(2)**4*y(3) + 1407*y(1)*y(2)**4*y(4) &
2999 & - 120*y(1)*y(2)**3*y(3)**2 - 120*y(1)*y(2)**3*y(3)*y(4) - 50*y(1)*y(2)**3*y(4) &
3000 & **2 + 150*y(1)*y(2)**2*y(3)**3 + 225*y(1)*y(2)**2*y(3)**2*y(4) + 125*y(1)*y(2) &
3001 & **2*y(3)*y(4)**2 + 25*y(1)*y(2)**2*y(4)**3 + 45*y(1)*y(2)*y(3)**4 + 90*y(1)*y(2) &
3002 & *y(3)**3*y(4) + 75*y(1)*y(2)*y(3)**2*y(4)**2 + 30*y(1)*y(2)*y(3)*y(4)**3 + 5*y(1) &
3003 & *y(2)*y(4)**4 + 3174*y(2)**6 + 6648*y(2)**5*y(3) + 3324*y(2)**5*y(4) + 4224*y(2) &
3004 & **4*y(3)**2 + 4224*y(2)**4*y(3)*y(4) + 1081*y(2)**4*y(4)**2 + 1200*y(2)**3*y(3) &
3005 & **3 + 1800*y(2)**3*y(3)**2*y(4) + 1000*y(2)**3*y(3)*y(4)**2 + 200*y(2)**3*y(4) &
3006 & **3 + 450*y(2)**2*y(3)**4 + 900*y(2)**2*y(3)**3*y(4) + 750*y(2)**2*y(3)**2*y(4) &
3007 & **2 + 300*y(2)**2*y(3)*y(4)**3 + 50*y(2)**2*y(4)**4))/(5*(y(2) + y(3))**2*(y(1) &
3008 & + y(2) + y(3))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
3010 & 4) = (4*y(2)**2*(105*y(1)**2*y(2)**3 + 220*y(1)**2*y(2)**2*y(3) + 110*y(1) &
3011 & **2*y(2)**2*y(4) + 35*y(1)**2*y(2)*y(3)**2 + 35*y(1)**2*y(2)*y(3)*y(4) + 5*y(1) &
3012 & **2*y(2)*y(4)**2 + 20*y(1)**2*y(3)**3 + 30*y(1)**2*y(3)**2*y(4) + 10*y(1)**2*y(3) &
3013 & *y(4)**2 - 1452*y(1)*y(2)**4 + 250*y(1)*y(2)**3*y(3) + 125*y(1)*y(2)**3*y(4) &
3014 & + 100*y(1)*y(2)**2*y(3)**2 + 100*y(1)*y(2)**2*y(3)*y(4) + 20*y(1)*y(2)**2*y(4) &
3015 & **2 + 90*y(1)*y(2)*y(3)**3 + 135*y(1)*y(2)*y(3)**2*y(4) + 55*y(1)*y(2)*y(3)*y(4) &
3016 & **2 + 5*y(1)*y(2)*y(4)**3 + 30*y(1)*y(3)**4 + 60*y(1)*y(3)**3*y(4) + 40*y(1)*y(3) &
3017 & **2*y(4)**2 + 10*y(1)*y(3)*y(4)**3 - 3219*y(2)**5 - 3694*y(2)**4*y(3) - 1847*y(2) &
3018 & **4*y(4) - 1040*y(2)**3*y(3)**2 - 1040*y(2)**3*y(3)*y(4) - 285*y(2)**3*y(4)**2 &
3019 & - 550*y(2)**2*y(3)**3 - 825*y(2)**2*y(3)**2*y(4) - 465*y(2)**2*y(3)*y(4)**2 &
3020 & - 95*y(2)**2*y(4)**3 + 15*y(2)*y(3)**4 + 30*y(2)*y(3)**3*y(4) + 20*y(2)*y(3) &
3021 & **2*y(4)**2 + 5*y(2)*y(3)*y(4)**3))/(5*(y(1) + y(2))*(y(2) + y(3))*(y(1) + y(2) &
3022 & + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
3024 & 5) = (4*y(2)**2*(831*y(2)**4 + 200*y(2)**3*y(3) + 100*y(2)**3*y(4) + 205*y(2) &
3025 & **2*y(3)**2 + 205*y(2)**2*y(3)*y(4) + 50*y(2)**2*y(4)**2 + 10*y(2)*y(3)**3 &
3026 & + 15*y(2)*y(3)**2*y(4) + 5*y(2)*y(3)*y(4)**2 + 5*y(3)**4 + 10*y(3)**3*y(4) &
3027 & + 5*y(3)**2*y(4)**2))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) &
3028 & + y(3) + y(4))**2)
3030 y = s_cb(i + 1:i + 4) - s_cb(i:i + 3)
3032 & 0) = (4*y(1)**2*(831*y(1)**4 + 200*y(1)**3*y(2) + 100*y(1)**3*y(3) + 205*y(1) &
3033 & **2*y(2)**2 + 205*y(1)**2*y(2)*y(3) + 50*y(1)**2*y(3)**2 + 10*y(1)*y(2)**3 &
3034 & + 15*y(1)*y(2)**2*y(3) + 5*y(1)*y(2)*y(3)**2 + 5*y(2)**4 + 10*y(2)**3*y(3) &
3035 & + 5*y(2)**2*y(3)**2))/(5*(y(3) + y(4))**2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) &
3036 & + y(3) + y(4))**2)
3038 & 1) = -(4*y(1)**2*(1662*y(1)**5 + 3824*y(1)**4*y(2) + 3624*y(1)**4*y(3) &
3039 & + 1762*y(1)**4*y(4) + 1515*y(1)**3*y(2)**2 + 2115*y(1)**3*y(2)*y(3) + 805*y(1) &
3040 & **3*y(2)*y(4) + 700*y(1)**3*y(3)**2 + 500*y(1)**3*y(3)*y(4) + 100*y(1)**3*y(4) &
3041 & **2 + 1060*y(1)**2*y(2)**3 + 2205*y(1)**2*y(2)**2*y(3) + 835*y(1)**2*y(2)**2*y(4) &
3042 & + 1445*y(1)**2*y(2)*y(3)**2 + 1030*y(1)**2*y(2)*y(3)*y(4) + 205*y(1)**2*y(2)*y(4) &
3043 & **2 + 300*y(1)**2*y(3)**3 + 300*y(1)**2*y(3)**2*y(4) + 100*y(1)**2*y(3)*y(4)**2 &
3044 & + 75*y(1)*y(2)**4 + 180*y(1)*y(2)**3*y(3) + 60*y(1)*y(2)**3*y(4) + 135*y(1)*y(2) &
3045 & **2*y(3)**2 + 90*y(1)*y(2)**2*y(3)*y(4) + 15*y(1)*y(2)**2*y(4)**2 + 30*y(1)*y(2) &
3046 & *y(3)**3 + 30*y(1)*y(2)*y(3)**2*y(4) + 10*y(1)*y(2)*y(3)*y(4)**2 + 30*y(2)**5 &
3047 & + 90*y(2)**4*y(3) + 30*y(2)**4*y(4) + 90*y(2)**3*y(3)**2 + 60*y(2)**3*y(3)*y(4) &
3048 & + 10*y(2)**3*y(4)**2 + 30*y(2)**2*y(3)**3 + 30*y(2)**2*y(3)**2*y(4) + 10*y(2) &
3049 & **2*y(3)*y(4)**2))/(5*(y(2) + y(3))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) &
3050 & + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
3052 & 2) = (4*y(1)**2*(1767*y(1)**4 + 725*y(1)**3*y(2) + 415*y(1)**3*y(3) + 105*y(4) &
3053 & *y(1)**3 + 665*y(1)**2*y(2)**2 + 775*y(1)**2*y(2)*y(3) + 220*y(4)*y(1)**2*y(2) &
3054 & + 215*y(1)**2*y(3)**2 + 110*y(4)*y(1)**2*y(3) + 75*y(1)*y(2)**3 + 130*y(1)*y(2) &
3055 & **2*y(3) + 35*y(4)*y(1)*y(2)**2 + 60*y(1)*y(2)*y(3)**2 + 35*y(4)*y(1)*y(2)*y(3) &
3056 & + 5*y(1)*y(3)**3 + 5*y(4)*y(1)*y(3)**2 + 30*y(2)**4 + 70*y(2)**3*y(3) + 20*y(4) &
3057 & *y(2)**3 + 50*y(2)**2*y(3)**2 + 30*y(4)*y(2)**2*y(3) + 10*y(2)*y(3)**3 + 10*y(4) &
3058 & *y(2)*y(3)**2))/(5*(y(1) + y(2))*(y(3) + y(4))*(y(1) + y(2) + y(3))*(y(2) + y(3) &
3059 & + y(4))*(y(1) + y(2) + y(3) + y(4))**2)
3061 & 3) = (4*y(1)**2*(831*y(1)**6 + 3624*y(1)**5*y(2) + 3524*y(1)**5*y(3) + 1762*y(1) &
3062 & **5*y(4) + 4884*y(1)**4*y(2)**2 + 9058*y(1)**4*y(2)*y(3) + 4529*y(1)**4*y(2)*y(4) &
3063 & + 4224*y(1)**4*y(3)**2 + 4224*y(1)**4*y(3)*y(4) + 1081*y(1)**4*y(4)**2 &
3064 & + 2565*y(1)**3*y(2)**3 + 6120*y(1)**3*y(2)**2*y(3) + 3060*y(1)**3*y(2)**2*y(4) &
3065 & + 4755*y(1)**3*y(2)*y(3)**2 + 4755*y(1)**3*y(2)*y(3)*y(4) + 1315*y(1)**3*y(2) &
3066 & *y(4)**2 + 1200*y(1)**3*y(3)**3 + 1800*y(1)**3*y(3)**2*y(4) + 1000*y(1)**3*y(3) &
3067 & *y(4)**2 + 200*y(1)**3*y(4)**3 + 1395*y(1)**2*y(2)**4 + 4380*y(1)**2*y(2)**3*y(3) &
3068 & + 2190*y(1)**2*y(2)**3*y(4) + 5025*y(1)**2*y(2)**2*y(3)**2 + 5025*y(1)**2*y(2) &
3069 & **2*y(3)*y(4) + 1390*y(1)**2*y(2)**2*y(4)**2 + 2490*y(1)**2*y(2)*y(3)**3 &
3070 & + 3735*y(1)**2*y(2)*y(3)**2*y(4) + 2075*y(1)**2*y(2)*y(3)*y(4)**2 + 415*y(1) &
3071 & **2*y(2)*y(4)**3 + 450*y(1)**2*y(3)**4 + 900*y(1)**2*y(3)**3*y(4) + 750*y(1) &
3072 & **2*y(3)**2*y(4)**2 + 300*y(1)**2*y(3)*y(4)**3 + 50*y(1)**2*y(4)**4 + 135*y(1) &
3073 & *y(2)**5 + 450*y(1)*y(2)**4*y(3) + 225*y(1)*y(2)**4*y(4) + 540*y(1)*y(2)**3*y(3) &
3074 & **2 + 540*y(1)*y(2)**3*y(3)*y(4) + 150*y(1)*y(2)**3*y(4)**2 + 270*y(1)*y(2) &
3075 & **2*y(3)**3 + 405*y(1)*y(2)**2*y(3)**2*y(4) + 225*y(1)*y(2)**2*y(3)*y(4)**2 &
3076 & + 45*y(1)*y(2)**2*y(4)**3 + 45*y(1)*y(2)*y(3)**4 + 90*y(1)*y(2)*y(3)**3*y(4) &
3077 & + 75*y(1)*y(2)*y(3)**2*y(4)**2 + 30*y(1)*y(2)*y(3)*y(4)**3 + 5*y(1)*y(2)*y(4)**4 &
3078 & + 45*y(2)**6 + 180*y(2)**5*y(3) + 90*y(2)**5*y(4) + 270*y(2)**4*y(3)**2 &
3079 & + 270*y(2)**4*y(3)*y(4) + 75*y(2)**4*y(4)**2 + 180*y(2)**3*y(3)**3 + 270*y(2) &
3080 & **3*y(3)**2*y(4) + 150*y(2)**3*y(3)*y(4)**2 + 30*y(2)**3*y(4)**3 + 45*y(2) &
3081 & **2*y(3)**4 + 90*y(2)**2*y(3)**3*y(4) + 75*y(2)**2*y(3)**2*y(4)**2 + 30*y(2) &
3082 & **2*y(3)*y(4)**3 + 5*y(2)**2*y(4)**4))/(5*(y(2) + y(3))**2*(y(1) + y(2) + y(3)) &
3083 & **2*(y(2) + y(3) + y(4))**2*(y(1) + y(2) + y(3) + y(4))**2)
3085 & 4) = -(4*y(1)**2*(1767*y(1)**5 + 4464*y(1)**4*y(2) + 4154*y(1)**4*y(3) &
3086 & + 2077*y(1)**4*y(4) + 2655*y(1)**3*y(2)**2 + 4010*y(1)**3*y(2)*y(3) + 2005*y(1) &
3087 & **3*y(2)*y(4) + 1460*y(1)**3*y(3)**2 + 1460*y(1)**3*y(3)*y(4) + 415*y(1)**3*y(4) &
3088 & **2 + 1800*y(1)**2*y(2)**3 + 4000*y(1)**2*y(2)**2*y(3) + 2000*y(1)**2*y(2) &
3089 & **2*y(4) + 2850*y(1)**2*y(2)*y(3)**2 + 2850*y(1)**2*y(2)*y(3)*y(4) + 790*y(1) &
3090 & **2*y(2)*y(4)**2 + 650*y(1)**2*y(3)**3 + 975*y(1)**2*y(3)**2*y(4) + 535*y(1) &
3091 & **2*y(3)*y(4)**2 + 105*y(1)**2*y(4)**3 + 270*y(1)*y(2)**4 + 720*y(1)*y(2)**3*y(3) &
3092 & + 360*y(1)*y(2)**3*y(4) + 645*y(1)*y(2)**2*y(3)**2 + 645*y(1)*y(2)**2*y(3)*y(4) &
3093 & + 165*y(1)*y(2)**2*y(4)**2 + 210*y(1)*y(2)*y(3)**3 + 315*y(1)*y(2)*y(3)**2*y(4) &
3094 & + 155*y(1)*y(2)*y(3)*y(4)**2 + 25*y(1)*y(2)*y(4)**3 + 15*y(1)*y(3)**4 + 30*y(1) &
3095 & *y(3)**3*y(4) + 20*y(1)*y(3)**2*y(4)**2 + 5*y(1)*y(3)*y(4)**3 + 90*y(2)**5 &
3096 & + 300*y(2)**4*y(3) + 150*y(2)**4*y(4) + 360*y(2)**3*y(3)**2 + 360*y(2)**3*y(3) &
3097 & *y(4) + 90*y(2)**3*y(4)**2 + 180*y(2)**2*y(3)**3 + 270*y(2)**2*y(3)**2*y(4) &
3098 & + 130*y(2)**2*y(3)*y(4)**2 + 20*y(2)**2*y(4)**3 + 30*y(2)*y(3)**4 + 60*y(2)*y(3) &
3099 & **3*y(4) + 40*y(2)*y(3)**2*y(4)**2 + 10*y(2)*y(3)*y(4)**3))/(5*(y(1) + y(2)) &
3100 & *(y(2) + y(3))*(y(1) + y(2) + y(3))**2*(y(2) + y(3) + y(4))*(y(1) + y(2) + y(3) &
3103 & 5) = (4*y(1)**2*(996*y(1)**4 + 675*y(1)**3*y(2) + 450*y(1)**3*y(3) + 225*y(1) &
3104 & **3*y(4) + 600*y(1)**2*y(2)**2 + 800*y(1)**2*y(2)*y(3) + 400*y(1)**2*y(2)*y(4) &
3105 & + 260*y(1)**2*y(3)**2 + 260*y(1)**2*y(3)*y(4) + 60*y(1)**2*y(4)**2 + 135*y(1) &
3106 & *y(2)**3 + 270*y(1)*y(2)**2*y(3) + 135*y(1)*y(2)**2*y(4) + 165*y(1)*y(2)*y(3)**2 &
3107 & + 165*y(1)*y(2)*y(3)*y(4) + 30*y(1)*y(2)*y(4)**2 + 30*y(1)*y(3)**3 + 45*y(1)*y(3) &
3108 & **2*y(4) + 15*y(1)*y(3)*y(4)**2 + 45*y(2)**4 + 120*y(2)**3*y(3) + 60*y(2)**3*y(4) &
3109 & + 110*y(2)**2*y(3)**2 + 110*y(2)**2*y(3)*y(4) + 20*y(2)**2*y(4)**2 + 40*y(2)*y(3) &
3110 & **3 + 60*y(2)*y(3)**2*y(4) + 20*y(2)*y(3)*y(4)**2 + 5*y(3)**4 + 10*y(3)**3*y(4) &
3111 & + 5*y(3)**2*y(4)**2))/(5*(y(1) + y(2))**2*(y(1) + y(2) + y(3))**2*(y(1) + y(2) &
3112 & + y(3) + y(4))**2)
3130# 866 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3136 h0 = (s_cb(s) - s_cb(0))/real(s, wp)
3138 if (abs((s_cb(i + 1) - s_cb(i)) - h0) > sqrt(epsilon(h0))*abs(h0))
then
3144 if (weno_dir == 1)
then
3146# 880 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3147#if defined(MFC_OpenACC)
3148# 880 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3150# 880 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3151#elif defined(MFC_OpenMP)
3152# 880 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3154# 880 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3156 else if (weno_dir == 2)
then
3158# 882 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3159#if defined(MFC_OpenACC)
3160# 882 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3162# 882 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3163#elif defined(MFC_OpenMP)
3164# 882 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3166# 882 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3170# 884 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3171#if defined(MFC_OpenACC)
3172# 884 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3174# 884 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3175#elif defined(MFC_OpenMP)
3176# 884 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"
3178# 884 "/home/runner/work/MFC/MFC/src/simulation/m_weno.fpp"