1% Vector test routine 2 3% Author: David Harper (algebra@liverpool.ac.uk) 4% Computer Algebra Support Officer 5% University of Liverpool Computer Laboratory. 6 7% Please compare carefully the output from running this test file with the 8% log file provided to make sure your implementation is correct. 9 10linelength 72; 11 12 1380 14 15off allfac; 16 17 on div; 18 19 20vec a,b,c; 21 22 23matrix q; 24 25 26a := avec(ax,ay,az); 27 28 29vec(x) := ax 30 31vec(y) := ay 32 33vec(z) := az 34 35 36b := avec(bx,by,bz); 37 38 39vec(x) := bx 40 41vec(y) := by 42 43vec(z) := bz 44 45 46q := mat((q11,q12,q13),(q21,q22,q23),(q31,q32,q33)); 47 48 49 [q11 q12 q13] 50 [ ] 51q := [q21 q22 q23] 52 [ ] 53 [q31 q32 q33] 54 55 56 57c := a+b; 58 59 60vec(x) := ax + bx 61 62vec(y) := ay + by 63 64vec(z) := az + bz 65 66 67c := a-b; 68 69 70vec(x) := ax - bx 71 72vec(y) := ay - by 73 74vec(z) := az - bz 75 76 77c := a cross b; 78 79 80vec(x) := ay*bz - az*by 81 82vec(y) := - ax*bz + az*bx 83 84vec(z) := ax*by - ay*bx 85 86 87d := a dot b; 88 89 90d := ax*bx + ay*by + az*bz 91 92a dot c; 93 94 950 96 b dot c; 97 98 990 100 101q*a; 102 103 104vec(x) := ax*q11 + ay*q12 + az*q13 105 106vec(y) := ax*q21 + ay*q22 + az*q23 107 108vec(z) := ax*q31 + ay*q32 + az*q33 109 110 111c:=2*f*a - b/7; 112 113 114 1 115vec(x) := 2*ax*f - ---*bx 116 7 117 118 1 119vec(y) := 2*ay*f - ---*by 120 7 121 122 1 123vec(z) := 2*az*f - ---*bz 124 7 125 126 127c(0); 128 129 130 1 1312*ax*f - ---*bx 132 7 133 c(1); 134 135 136 1 1372*ay*f - ---*by 138 7 139 c(2); 140 141 142 1 1432*az*f - ---*bz 144 7 145 146 1471/vmod(a); 148 149 150 2 2 2 - 1/2 151(ax + ay + az ) 152 153b/vmod(a); 154 155 156 2 2 2 - 1/2 157vec(x) := (ax + ay + az ) *bx 158 159 2 2 2 - 1/2 160vec(y) := (ax + ay + az ) *by 161 162 2 2 2 - 1/2 163vec(z) := (ax + ay + az ) *bz 164 165 166(a cross b)/(a dot b); 167 168 169 ay*bz - az*by 170vec(x) := ----------------------- 171 ax*bx + ay*by + az*bz 172 173 - ax*bz + az*bx 174vec(y) := ----------------------- 175 ax*bx + ay*by + az*bz 176 177 ax*by - ay*bx 178vec(z) := ----------------------- 179 ax*bx + ay*by + az*bz 180 181 1822/3*vmod(a)*a*(a dot c)/(vmod(a cross c)); 183 184 185 28 2 2 2 2 186vec(x) := ----*(ax *by + ax *bz - 2*ax*ay*bx*by - 2*ax*az*bx*bz 187 3 188 189 2 2 2 2 2 2 2 2 190 + ay *bx + ay *bz - 2*ay*az*by*bz + az *bx + az *by 191 192 - 1 2 2 2 3 2 2 2 193 )**------*sqrt(ax + ay + az )*ax *f - ---*(ax *by 194 2 3 195 196 2 2 2 2 197 + ax *bz - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx 198 199 2 2 2 2 2 2 - 1 200 + ay *bz - 2*ay*az*by*bz + az *bx + az *by )**------ 201 2 202 203 2 2 2 2 28 2 2 2 2 204 *sqrt(ax + ay + az )*ax *bx + ----*(ax *by + ax *bz 205 3 206 207 2 2 2 2 208 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 209 210 2 2 2 2 - 1 211 - 2*ay*az*by*bz + az *bx + az *by )**------ 212 2 213 214 2 2 2 2 2 2 2 2 2 215 *sqrt(ax + ay + az )*ax*ay *f - ---*(ax *by + ax *bz 216 3 217 218 2 2 2 2 219 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 220 221 2 2 2 2 - 1 222 - 2*ay*az*by*bz + az *bx + az *by )**------ 223 2 224 225 2 2 2 28 2 2 2 2 226 *sqrt(ax + ay + az )*ax*ay*by + ----*(ax *by + ax *bz 227 3 228 229 2 2 2 2 230 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 231 232 2 2 2 2 - 1 233 - 2*ay*az*by*bz + az *bx + az *by )**------ 234 2 235 236 2 2 2 2 2 2 2 2 2 237 *sqrt(ax + ay + az )*ax*az *f - ---*(ax *by + ax *bz 238 3 239 240 2 2 2 2 241 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 242 243 2 2 2 2 - 1 244 - 2*ay*az*by*bz + az *bx + az *by )**------ 245 2 246 247 2 2 2 248 *sqrt(ax + ay + az )*ax*az*bz 249 250 28 2 2 2 2 251vec(y) := ----*(ax *by + ax *bz - 2*ax*ay*bx*by - 2*ax*az*bx*bz 252 3 253 254 2 2 2 2 2 2 2 2 255 + ay *bx + ay *bz - 2*ay*az*by*bz + az *bx + az *by 256 257 - 1 2 2 2 2 2 2 2 258 )**------*sqrt(ax + ay + az )*ax *ay*f - ---*(ax *by 259 2 3 260 261 2 2 2 2 262 + ax *bz - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx 263 264 2 2 2 2 2 2 - 1 265 + ay *bz - 2*ay*az*by*bz + az *bx + az *by )**------ 266 2 267 268 2 2 2 28 2 2 2 2 269 *sqrt(ax + ay + az )*ax*ay*bx + ----*(ax *by + ax *bz 270 3 271 272 2 2 2 2 273 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 274 275 2 2 2 2 - 1 276 - 2*ay*az*by*bz + az *bx + az *by )**------ 277 2 278 279 2 2 2 3 2 2 2 2 2 280 *sqrt(ax + ay + az )*ay *f - ---*(ax *by + ax *bz 281 3 282 283 2 2 2 2 284 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 285 286 2 2 2 2 - 1 287 - 2*ay*az*by*bz + az *bx + az *by )**------ 288 2 289 290 2 2 2 2 28 2 2 2 2 291 *sqrt(ax + ay + az )*ay *by + ----*(ax *by + ax *bz 292 3 293 294 2 2 2 2 295 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 296 297 2 2 2 2 - 1 298 - 2*ay*az*by*bz + az *bx + az *by )**------ 299 2 300 301 2 2 2 2 2 2 2 2 2 302 *sqrt(ax + ay + az )*ay*az *f - ---*(ax *by + ax *bz 303 3 304 305 2 2 2 2 306 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 307 308 2 2 2 2 - 1 309 - 2*ay*az*by*bz + az *bx + az *by )**------ 310 2 311 312 2 2 2 313 *sqrt(ax + ay + az )*ay*az*bz 314 315 28 2 2 2 2 316vec(z) := ----*(ax *by + ax *bz - 2*ax*ay*bx*by - 2*ax*az*bx*bz 317 3 318 319 2 2 2 2 2 2 2 2 320 + ay *bx + ay *bz - 2*ay*az*by*bz + az *bx + az *by 321 322 - 1 2 2 2 2 2 2 2 323 )**------*sqrt(ax + ay + az )*ax *az*f - ---*(ax *by 324 2 3 325 326 2 2 2 2 327 + ax *bz - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx 328 329 2 2 2 2 2 2 - 1 330 + ay *bz - 2*ay*az*by*bz + az *bx + az *by )**------ 331 2 332 333 2 2 2 28 2 2 2 2 334 *sqrt(ax + ay + az )*ax*az*bx + ----*(ax *by + ax *bz 335 3 336 337 2 2 2 2 338 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 339 340 2 2 2 2 - 1 341 - 2*ay*az*by*bz + az *bx + az *by )**------ 342 2 343 344 2 2 2 2 2 2 2 2 2 345 *sqrt(ax + ay + az )*ay *az*f - ---*(ax *by + ax *bz 346 3 347 348 2 2 2 2 349 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 350 351 2 2 2 2 - 1 352 - 2*ay*az*by*bz + az *bx + az *by )**------ 353 2 354 355 2 2 2 28 2 2 2 2 356 *sqrt(ax + ay + az )*ay*az*by + ----*(ax *by + ax *bz 357 3 358 359 2 2 2 2 360 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 361 362 2 2 2 2 - 1 363 - 2*ay*az*by*bz + az *bx + az *by )**------ 364 2 365 366 2 2 2 3 2 2 2 2 2 367 *sqrt(ax + ay + az )*az *f - ---*(ax *by + ax *bz 368 3 369 370 2 2 2 2 371 - 2*ax*ay*bx*by - 2*ax*az*bx*bz + ay *bx + ay *bz 372 373 2 2 2 2 - 1 374 - 2*ay*az*by*bz + az *bx + az *by )**------ 375 2 376 377 2 2 2 2 378 *sqrt(ax + ay + az )*az *bz 379 380 381 382a := avec(x**2*y**3,log(z+x),13*z-y); 383 384 385 2 3 386vec(x) := x *y 387 388vec(y) := log(x + z) 389 390vec(z) := - y + 13*z 391 392 393df(a,x); 394 395 396 3 397vec(x) := 2*x*y 398 399 1 400vec(y) := ------- 401 x + z 402 403vec(z) := 0 404 405 406df(a,x,y); 407 408 409 2 410vec(x) := 6*x*y 411 412vec(y) := 0 413 414vec(z) := 0 415 416 417int(a,x); 418 419 420 1 3 3 421vec(x) := ---*x *y 422 3 423 424vec(y) := log(x + z)*x + log(x + z)*z - x 425 426vec(z) := - x*y + 13*x*z 427 428 429exp(a); 430 431 432 2 3 433 x *y 434vec(x) := e 435 436vec(y) := x + z 437 438 - y + 13*z 439vec(z) := e 440 441 442log sin b; 443 444 445vec(x) := log(sin(bx)) 446 447vec(y) := log(sin(by)) 448 449vec(z) := log(sin(bz)) 450 451 452 453a := avec(ax,ay,az); 454 455 456vec(x) := ax 457 458vec(y) := ay 459 460vec(z) := az 461 462 463depend ax,x,y,z; 464 465 depend ay,x,y,z; 466 467 depend az,x,y,z; 468 469 470depend p,x,y,z; 471 472 473c := grad p; 474 475 476vec(x) := df(p,x) 477 478vec(y) := df(p,y) 479 480vec(z) := df(p,z) 481 482 483div c; 484 485 486df(p,x,2) + df(p,y,2) + df(p,z,2) 487 488delsq p; 489 490 491df(p,x,2) + df(p,y,2) + df(p,z,2) 492 493div a; 494 495 496df(ax,x) + df(ay,y) + df(az,z) 497 498curl a; 499 500 501vec(x) := - df(ay,z) + df(az,y) 502 503vec(y) := df(ax,z) - df(az,x) 504 505vec(z) := - df(ax,y) + df(ay,x) 506 507 508delsq a; 509 510 511vec(x) := df(ax,x,2) + df(ax,y,2) + df(ax,z,2) 512 513vec(y) := df(ay,x,2) + df(ay,y,2) + df(ay,z,2) 514 515vec(z) := df(az,x,2) + df(az,y,2) + df(az,z,2) 516 517 518 519depend h1,x,y,z; 520 521 depend h2,x,y,z; 522 523 depend h3,x,y,z; 524 525 526scalefactors(h1,h2,h3); 527 528 529grad p; 530 531 532 -1 533vec(x) := df(p,x)*h1 534 535 -1 536vec(y) := df(p,y)*h2 537 538 -1 539vec(z) := df(p,z)*h3 540 541 542div a; 543 544 545 -1 -1 -1 -1 -1 546df(ax,x)*h1 + df(ay,y)*h2 + df(az,z)*h3 + df(h1,y)*ay*h1 *h2 547 548 -1 -1 -1 -1 -1 -1 549 + df(h1,z)*az*h1 *h3 + df(h2,x)*ax*h1 *h2 + df(h2,z)*az*h2 *h3 550 551 -1 -1 -1 -1 552 + df(h3,x)*ax*h1 *h3 + df(h3,y)*ay*h2 *h3 553 554curl a; 555 556 557 -1 -1 -1 -1 558vec(x) := - df(ay,z)*h3 + df(az,y)*h2 - df(h2,z)*ay*h2 *h3 559 560 -1 -1 561 + df(h3,y)*az*h2 *h3 562 563 -1 -1 -1 -1 564vec(y) := df(ax,z)*h3 - df(az,x)*h1 + df(h1,z)*ax*h1 *h3 565 566 -1 -1 567 - df(h3,x)*az*h1 *h3 568 569 -1 -1 -1 -1 570vec(z) := - df(ax,y)*h2 + df(ay,x)*h1 - df(h1,y)*ax*h1 *h2 571 572 -1 -1 573 + df(h2,x)*ay*h1 *h2 574 575 576dp1 := delsq p; 577 578 579 -3 -1 -2 580dp1 := - df(h1,x)*df(p,x)*h1 + df(h1,y)*df(p,y)*h1 *h2 581 582 -1 -2 -2 -1 583 + df(h1,z)*df(p,z)*h1 *h3 + df(h2,x)*df(p,x)*h1 *h2 584 585 -3 -1 -2 586 - df(h2,y)*df(p,y)*h2 + df(h2,z)*df(p,z)*h2 *h3 587 588 -2 -1 -2 -1 589 + df(h3,x)*df(p,x)*h1 *h3 + df(h3,y)*df(p,y)*h2 *h3 590 591 -3 -2 -2 592 - df(h3,z)*df(p,z)*h3 + df(p,x,2)*h1 + df(p,y,2)*h2 593 594 -2 595 + df(p,z,2)*h3 596 597dp2 := div grad p; 598 599 600 -3 -1 -2 601dp2 := - df(h1,x)*df(p,x)*h1 + df(h1,y)*df(p,y)*h1 *h2 602 603 -1 -2 -2 -1 604 + df(h1,z)*df(p,z)*h1 *h3 + df(h2,x)*df(p,x)*h1 *h2 605 606 -3 -1 -2 607 - df(h2,y)*df(p,y)*h2 + df(h2,z)*df(p,z)*h2 *h3 608 609 -2 -1 -2 -1 610 + df(h3,x)*df(p,x)*h1 *h3 + df(h3,y)*df(p,y)*h2 *h3 611 612 -3 -2 -2 613 - df(h3,z)*df(p,z)*h3 + df(p,x,2)*h1 + df(p,y,2)*h2 614 615 -2 616 + df(p,z,2)*h3 617 618dp1-dp2; 619 620 6210 622 623delsq a; 624 625 626 -2 -3 627vec(x) := df(ax,x,2)*h1 - df(ax,x)*df(h1,x)*h1 628 629 -2 -1 -2 -1 630 + df(ax,x)*df(h2,x)*h1 *h2 + df(ax,x)*df(h3,x)*h1 *h3 631 632 -2 -1 -2 633 + df(ax,y,2)*h2 + df(ax,y)*df(h1,y)*h1 *h2 634 635 -3 -2 -1 636 - df(ax,y)*df(h2,y)*h2 + df(ax,y)*df(h3,y)*h2 *h3 637 638 -2 -1 -2 639 + df(ax,z,2)*h3 + df(ax,z)*df(h1,z)*h1 *h3 640 641 -1 -2 -3 642 + df(ax,z)*df(h2,z)*h2 *h3 - df(ax,z)*df(h3,z)*h3 643 644 -2 -1 645 + 2*df(ay,x)*df(h1,y)*h1 *h2 646 647 -1 -2 648 - 2*df(ay,y)*df(h2,x)*h1 *h2 649 650 -2 -1 651 + 2*df(az,x)*df(h1,z)*h1 *h3 652 653 -1 -2 -2 -1 654 - 2*df(az,z)*df(h3,x)*h1 *h3 + df(h1,x,y)*ay*h1 *h2 655 656 -2 -1 -3 -1 657 + df(h1,x,z)*az*h1 *h3 - df(h1,x)*df(h1,y)*ay*h1 *h2 658 659 -3 -1 660 - df(h1,x)*df(h1,z)*az*h1 *h3 661 662 -3 -1 663 - df(h1,x)*df(h2,x)*ax*h1 *h2 664 665 -3 -1 -1 -2 666 - df(h1,x)*df(h3,x)*ax*h1 *h3 + df(h1,y,2)*ax*h1 *h2 667 668 2 -2 -2 -1 -3 669 - df(h1,y) *ax*h1 *h2 - df(h1,y)*df(h2,y)*ax*h1 *h2 670 671 -1 -2 -1 672 + df(h1,y)*df(h3,y)*ax*h1 *h2 *h3 673 674 -1 -2 2 -2 -2 675 + df(h1,z,2)*ax*h1 *h3 - df(h1,z) *ax*h1 *h3 676 677 -1 -1 -2 678 + df(h1,z)*df(h2,z)*ax*h1 *h2 *h3 679 680 -1 -3 -1 -2 681 - df(h1,z)*df(h3,z)*ax*h1 *h3 - df(h2,x,y)*ay*h1 *h2 682 683 -1 -1 -1 -2 -1 684 + df(h2,x,z)*az*h1 *h2 *h3 + df(h2,x,2)*ax*h1 *h2 685 686 2 -2 -2 -1 -3 687 - df(h2,x) *ax*h1 *h2 + df(h2,x)*df(h2,y)*ay*h1 *h2 688 689 -1 -2 -1 690 - df(h2,x)*df(h2,z)*az*h1 *h2 *h3 691 692 -1 -2 -1 693 - 2*df(h2,x)*df(h3,y)*ay*h1 *h2 *h3 694 695 -1 -1 -2 696 - 2*df(h2,z)*df(h3,x)*az*h1 *h2 *h3 697 698 -1 -1 -1 -1 -2 699 + df(h3,x,y)*ay*h1 *h2 *h3 - df(h3,x,z)*az*h1 *h3 700 701 -2 -1 2 -2 -2 702 + df(h3,x,2)*ax*h1 *h3 - df(h3,x) *ax*h1 *h3 703 704 -1 -1 -2 705 - df(h3,x)*df(h3,y)*ay*h1 *h2 *h3 706 707 -1 -3 708 + df(h3,x)*df(h3,z)*az*h1 *h3 709 710 -2 -1 711vec(y) := - 2*df(ax,x)*df(h1,y)*h1 *h2 712 713 -1 -2 -2 714 + 2*df(ax,y)*df(h2,x)*h1 *h2 + df(ay,x,2)*h1 715 716 -3 -2 -1 717 - df(ay,x)*df(h1,x)*h1 + df(ay,x)*df(h2,x)*h1 *h2 718 719 -2 -1 -2 720 + df(ay,x)*df(h3,x)*h1 *h3 + df(ay,y,2)*h2 721 722 -1 -2 -3 723 + df(ay,y)*df(h1,y)*h1 *h2 - df(ay,y)*df(h2,y)*h2 724 725 -2 -1 -2 726 + df(ay,y)*df(h3,y)*h2 *h3 + df(ay,z,2)*h3 727 728 -1 -2 -1 -2 729 + df(ay,z)*df(h1,z)*h1 *h3 + df(ay,z)*df(h2,z)*h2 *h3 730 731 -3 -2 -1 732 - df(ay,z)*df(h3,z)*h3 + 2*df(az,y)*df(h2,z)*h2 *h3 733 734 -1 -2 -2 -1 735 - 2*df(az,z)*df(h3,y)*h2 *h3 - df(h1,x,y)*ax*h1 *h2 736 737 -3 -1 738 + df(h1,x)*df(h1,y)*ax*h1 *h2 739 740 -3 -1 741 - df(h1,x)*df(h2,x)*ay*h1 *h2 742 743 -1 -1 -1 -1 -2 744 + df(h1,y,z)*az*h1 *h2 *h3 + df(h1,y,2)*ay*h1 *h2 745 746 2 -2 -2 747 - df(h1,y) *ay*h1 *h2 748 749 -2 -1 -1 750 - df(h1,y)*df(h1,z)*az*h1 *h2 *h3 751 752 -1 -3 753 - df(h1,y)*df(h2,y)*ay*h1 *h2 754 755 -2 -1 -1 756 - 2*df(h1,y)*df(h3,x)*ax*h1 *h2 *h3 757 758 -1 -1 -2 759 + df(h1,z)*df(h2,z)*ay*h1 *h2 *h3 760 761 -1 -1 -2 762 - 2*df(h1,z)*df(h3,y)*az*h1 *h2 *h3 763 764 -1 -2 -2 -1 765 + df(h2,x,y)*ax*h1 *h2 + df(h2,x,2)*ay*h1 *h2 766 767 2 -2 -2 -1 -3 768 - df(h2,x) *ay*h1 *h2 - df(h2,x)*df(h2,y)*ax*h1 *h2 769 770 -2 -1 -1 771 + df(h2,x)*df(h3,x)*ay*h1 *h2 *h3 772 773 -2 -1 -3 -1 774 + df(h2,y,z)*az*h2 *h3 - df(h2,y)*df(h2,z)*az*h2 *h3 775 776 -3 -1 -1 -2 777 - df(h2,y)*df(h3,y)*ay*h2 *h3 + df(h2,z,2)*ay*h2 *h3 778 779 2 -2 -2 -1 -3 780 - df(h2,z) *ay*h2 *h3 - df(h2,z)*df(h3,z)*ay*h2 *h3 781 782 -1 -1 -1 783 + df(h3,x,y)*ax*h1 *h2 *h3 784 785 -1 -1 -2 786 - df(h3,x)*df(h3,y)*ax*h1 *h2 *h3 787 788 -1 -2 -2 -1 789 - df(h3,y,z)*az*h2 *h3 + df(h3,y,2)*ay*h2 *h3 790 791 2 -2 -2 -1 -3 792 - df(h3,y) *ay*h2 *h3 + df(h3,y)*df(h3,z)*az*h2 *h3 793 794 -2 -1 795vec(z) := - 2*df(ax,x)*df(h1,z)*h1 *h3 796 797 -1 -2 798 + 2*df(ax,z)*df(h3,x)*h1 *h3 799 800 -2 -1 801 - 2*df(ay,y)*df(h2,z)*h2 *h3 802 803 -1 -2 -2 804 + 2*df(ay,z)*df(h3,y)*h2 *h3 + df(az,x,2)*h1 805 806 -3 -2 -1 807 - df(az,x)*df(h1,x)*h1 + df(az,x)*df(h2,x)*h1 *h2 808 809 -2 -1 -2 810 + df(az,x)*df(h3,x)*h1 *h3 + df(az,y,2)*h2 811 812 -1 -2 -3 813 + df(az,y)*df(h1,y)*h1 *h2 - df(az,y)*df(h2,y)*h2 814 815 -2 -1 -2 816 + df(az,y)*df(h3,y)*h2 *h3 + df(az,z,2)*h3 817 818 -1 -2 -1 -2 819 + df(az,z)*df(h1,z)*h1 *h3 + df(az,z)*df(h2,z)*h2 *h3 820 821 -3 -2 -1 822 - df(az,z)*df(h3,z)*h3 - df(h1,x,z)*ax*h1 *h3 823 824 -3 -1 825 + df(h1,x)*df(h1,z)*ax*h1 *h3 826 827 -3 -1 828 - df(h1,x)*df(h3,x)*az*h1 *h3 829 830 -1 -1 -1 831 + df(h1,y,z)*ay*h1 *h2 *h3 832 833 -2 -1 -1 834 - df(h1,y)*df(h1,z)*ay*h1 *h2 *h3 835 836 -1 -2 -1 837 - 2*df(h1,y)*df(h2,z)*ay*h1 *h2 *h3 838 839 -1 -2 -1 840 + df(h1,y)*df(h3,y)*az*h1 *h2 *h3 841 842 -1 -2 2 -2 -2 843 + df(h1,z,2)*az*h1 *h3 - df(h1,z) *az*h1 *h3 844 845 -2 -1 -1 846 - 2*df(h1,z)*df(h2,x)*ax*h1 *h2 *h3 847 848 -1 -3 849 - df(h1,z)*df(h3,z)*az*h1 *h3 850 851 -1 -1 -1 852 + df(h2,x,z)*ax*h1 *h2 *h3 853 854 -1 -2 -1 855 - df(h2,x)*df(h2,z)*ax*h1 *h2 *h3 856 857 -2 -1 -1 858 + df(h2,x)*df(h3,x)*az*h1 *h2 *h3 859 860 -2 -1 -3 -1 861 - df(h2,y,z)*ay*h2 *h3 + df(h2,y)*df(h2,z)*ay*h2 *h3 862 863 -3 -1 -1 -2 864 - df(h2,y)*df(h3,y)*az*h2 *h3 + df(h2,z,2)*az*h2 *h3 865 866 2 -2 -2 -1 -3 867 - df(h2,z) *az*h2 *h3 - df(h2,z)*df(h3,z)*az*h2 *h3 868 869 -1 -2 -2 -1 870 + df(h3,x,z)*ax*h1 *h3 + df(h3,x,2)*az*h1 *h3 871 872 2 -2 -2 -1 -3 873 - df(h3,x) *az*h1 *h3 - df(h3,x)*df(h3,z)*ax*h1 *h3 874 875 -1 -2 -2 -1 876 + df(h3,y,z)*ay*h2 *h3 + df(h3,y,2)*az*h2 *h3 877 878 2 -2 -2 -1 -3 879 - df(h3,y) *az*h2 *h3 - df(h3,y)*df(h3,z)*ay*h2 *h3 880 881 882curl grad p; 883 884 885vec(x) := 0 886 887vec(y) := 0 888 889vec(z) := 0 890 891 892grad div a; 893 894 895 -2 -3 896vec(x) := df(ax,x,2)*h1 - df(ax,x)*df(h1,x)*h1 897 898 -2 -1 -2 -1 899 + df(ax,x)*df(h2,x)*h1 *h2 + df(ax,x)*df(h3,x)*h1 *h3 900 901 -1 -1 -2 -1 902 + df(ay,x,y)*h1 *h2 + df(ay,x)*df(h1,y)*h1 *h2 903 904 -1 -1 -1 905 + df(ay,x)*df(h3,y)*h1 *h2 *h3 906 907 -1 -2 -1 -1 908 - df(ay,y)*df(h2,x)*h1 *h2 + df(az,x,z)*h1 *h3 909 910 -2 -1 911 + df(az,x)*df(h1,z)*h1 *h3 912 913 -1 -1 -1 914 + df(az,x)*df(h2,z)*h1 *h2 *h3 915 916 -1 -2 -2 -1 917 - df(az,z)*df(h3,x)*h1 *h3 + df(h1,x,y)*ay*h1 *h2 918 919 -2 -1 -3 -1 920 + df(h1,x,z)*az*h1 *h3 - df(h1,x)*df(h1,y)*ay*h1 *h2 921 922 -3 -1 923 - df(h1,x)*df(h1,z)*az*h1 *h3 924 925 -3 -1 926 - df(h1,x)*df(h2,x)*ax*h1 *h2 927 928 -3 -1 929 - df(h1,x)*df(h3,x)*ax*h1 *h3 930 931 -2 -2 932 - df(h1,y)*df(h2,x)*ay*h1 *h2 933 934 -2 -2 935 - df(h1,z)*df(h3,x)*az*h1 *h3 936 937 -1 -1 -1 -2 -1 938 + df(h2,x,z)*az*h1 *h2 *h3 + df(h2,x,2)*ax*h1 *h2 939 940 2 -2 -2 941 - df(h2,x) *ax*h1 *h2 942 943 -1 -2 -1 944 - df(h2,x)*df(h2,z)*az*h1 *h2 *h3 945 946 -1 -2 -1 947 - df(h2,x)*df(h3,y)*ay*h1 *h2 *h3 948 949 -1 -1 -2 950 - df(h2,z)*df(h3,x)*az*h1 *h2 *h3 951 952 -1 -1 -1 -2 -1 953 + df(h3,x,y)*ay*h1 *h2 *h3 + df(h3,x,2)*ax*h1 *h3 954 955 2 -2 -2 956 - df(h3,x) *ax*h1 *h3 957 958 -1 -1 -2 959 - df(h3,x)*df(h3,y)*ay*h1 *h2 *h3 960 961 -1 -1 -2 -1 962vec(y) := df(ax,x,y)*h1 *h2 - df(ax,x)*df(h1,y)*h1 *h2 963 964 -1 -2 965 + df(ax,y)*df(h2,x)*h1 *h2 966 967 -1 -1 -1 -2 968 + df(ax,y)*df(h3,x)*h1 *h2 *h3 + df(ay,y,2)*h2 969 970 -1 -2 -3 971 + df(ay,y)*df(h1,y)*h1 *h2 - df(ay,y)*df(h2,y)*h2 972 973 -2 -1 -1 -1 974 + df(ay,y)*df(h3,y)*h2 *h3 + df(az,y,z)*h2 *h3 975 976 -1 -1 -1 977 + df(az,y)*df(h1,z)*h1 *h2 *h3 978 979 -2 -1 -1 -2 980 + df(az,y)*df(h2,z)*h2 *h3 - df(az,z)*df(h3,y)*h2 *h3 981 982 -1 -1 -1 -1 -2 983 + df(h1,y,z)*az*h1 *h2 *h3 + df(h1,y,2)*ay*h1 *h2 984 985 2 -2 -2 986 - df(h1,y) *ay*h1 *h2 987 988 -2 -1 -1 989 - df(h1,y)*df(h1,z)*az*h1 *h2 *h3 990 991 -2 -2 992 - df(h1,y)*df(h2,x)*ax*h1 *h2 993 994 -1 -3 995 - df(h1,y)*df(h2,y)*ay*h1 *h2 996 997 -2 -1 -1 998 - df(h1,y)*df(h3,x)*ax*h1 *h2 *h3 999 1000 -1 -1 -2 1001 - df(h1,z)*df(h3,y)*az*h1 *h2 *h3 1002 1003 -1 -2 -1 -3 1004 + df(h2,x,y)*ax*h1 *h2 - df(h2,x)*df(h2,y)*ax*h1 *h2 1005 1006 -2 -1 -3 -1 1007 + df(h2,y,z)*az*h2 *h3 - df(h2,y)*df(h2,z)*az*h2 *h3 1008 1009 -3 -1 1010 - df(h2,y)*df(h3,y)*ay*h2 *h3 1011 1012 -2 -2 1013 - df(h2,z)*df(h3,y)*az*h2 *h3 1014 1015 -1 -1 -1 1016 + df(h3,x,y)*ax*h1 *h2 *h3 1017 1018 -1 -1 -2 1019 - df(h3,x)*df(h3,y)*ax*h1 *h2 *h3 1020 1021 -2 -1 2 -2 -2 1022 + df(h3,y,2)*ay*h2 *h3 - df(h3,y) *ay*h2 *h3 1023 1024 -1 -1 -2 -1 1025vec(z) := df(ax,x,z)*h1 *h3 - df(ax,x)*df(h1,z)*h1 *h3 1026 1027 -1 -1 -1 1028 + df(ax,z)*df(h2,x)*h1 *h2 *h3 1029 1030 -1 -2 -1 -1 1031 + df(ax,z)*df(h3,x)*h1 *h3 + df(ay,y,z)*h2 *h3 1032 1033 -2 -1 1034 - df(ay,y)*df(h2,z)*h2 *h3 1035 1036 -1 -1 -1 1037 + df(ay,z)*df(h1,y)*h1 *h2 *h3 1038 1039 -1 -2 -2 1040 + df(ay,z)*df(h3,y)*h2 *h3 + df(az,z,2)*h3 1041 1042 -1 -2 -1 -2 1043 + df(az,z)*df(h1,z)*h1 *h3 + df(az,z)*df(h2,z)*h2 *h3 1044 1045 -3 -1 -1 -1 1046 - df(az,z)*df(h3,z)*h3 + df(h1,y,z)*ay*h1 *h2 *h3 1047 1048 -2 -1 -1 1049 - df(h1,y)*df(h1,z)*ay*h1 *h2 *h3 1050 1051 -1 -2 -1 1052 - df(h1,y)*df(h2,z)*ay*h1 *h2 *h3 1053 1054 -1 -2 2 -2 -2 1055 + df(h1,z,2)*az*h1 *h3 - df(h1,z) *az*h1 *h3 1056 1057 -2 -1 -1 1058 - df(h1,z)*df(h2,x)*ax*h1 *h2 *h3 1059 1060 -2 -2 1061 - df(h1,z)*df(h3,x)*ax*h1 *h3 1062 1063 -1 -3 1064 - df(h1,z)*df(h3,z)*az*h1 *h3 1065 1066 -1 -1 -1 1067 + df(h2,x,z)*ax*h1 *h2 *h3 1068 1069 -1 -2 -1 1070 - df(h2,x)*df(h2,z)*ax*h1 *h2 *h3 1071 1072 -1 -2 2 -2 -2 1073 + df(h2,z,2)*az*h2 *h3 - df(h2,z) *az*h2 *h3 1074 1075 -2 -2 1076 - df(h2,z)*df(h3,y)*ay*h2 *h3 1077 1078 -1 -3 -1 -2 1079 - df(h2,z)*df(h3,z)*az*h2 *h3 + df(h3,x,z)*ax*h1 *h3 1080 1081 -1 -3 -1 -2 1082 - df(h3,x)*df(h3,z)*ax*h1 *h3 + df(h3,y,z)*ay*h2 *h3 1083 1084 -1 -3 1085 - df(h3,y)*df(h3,z)*ay*h2 *h3 1086 1087 1088div curl a; 1089 1090 10910 1092 1093 1094% Examples of integration : (1) Volume integrals 1095 1096getcsystem 'spherical; 1097 1098 1099(r theta phi) 1100 1101 1102% Example 1 : integration of r**n over a sphere 1103 1104origin := avec(0,0,0); 1105 1106 1107vec(r) := 0 1108 1109vec(theta) := 0 1110 1111vec(phi) := 0 1112 1113 1114upperbound := avec(rr,pi,2*pi); 1115 1116 1117vec(r) := rr 1118 1119vec(theta) := pi 1120 1121vec(phi) := 2*pi 1122 1123 1124volintegral(r**n,origin,upperbound); 1125 1126 1127 n 3 1128 4*rr *pi*rr 1129-------------- 1130 n + 3 1131 1132 1133% Substitute n=0 to get the volume of a sphere 1134sub(n=0,ws); 1135 1136 1137 4 3 1138---*pi*rr 1139 3 1140 1141 1142% Example 2 : volume of a right-circular cone 1143 1144getcsystem 'cylindrical; 1145 1146 1147(r z phi) 1148 1149upperbound := avec(pp*z,h,2*pi); 1150 1151 1152vec(r) := pp*z 1153 1154vec(z) := h 1155 1156vec(phi) := 2*pi 1157 1158 1159volintorder := avec(2,0,1); 1160 1161 1162vec(r) := 2 1163 1164vec(z) := 0 1165 1166vec(phi) := 1 1167 1168 % Integrate in the order : phi, r, z 1169cone := volintegral(1,origin,upperbound); 1170 1171 1172 1 3 2 1173cone := ---*h *pi*pp 1174 3 1175 1176 1177% Now we replace P*Z by RR to get the result in the familiar form 1178 1179let pp*h=rr; 1180 1181 1182cone := cone; 1183 1184 1185 1 2 1186cone := ---*h*pi*rr 1187 3 1188 % This is the familiar form 1189clear pp*h; 1190 1191 1192 1193% Example 3 : line integral to obtain the length of a line of latitude 1194% on a sphere 1195 1196getcsystem 'spherical; 1197 1198 1199(r theta phi) 1200 1201 1202a := avec(0,0,1); 1203 1204 1205vec(r) := 0 1206 1207vec(theta) := 0 1208 1209vec(phi) := 1 1210 1211 % Function vector is the tangent to the 1212 % line of latitude 1213curve := avec(rr,latitude,phi); 1214 1215 1216vec(r) := rr 1217 1218vec(theta) := latitude 1219 1220vec(phi) := phi 1221 1222 % Path is round a line of latitude 1223 1224deflineint(a,curve,phi,0,2*pi); 1225 1226 12272*sin(latitude)*pi*rr 1228 1229 1230end; 1231 1232Tested on x86_64-pc-windows CSL 1233Time (counter 1): 0 ms plus GC time: 16 ms 1234 1235End of Lisp run after 0.00+0.06 seconds 1236real 0.22 1237user 0.01 1238sys 0.04 1239