%  
%   pfunction.m   ver 1.0  October 16, 2005  
%
%   by Tom Irvine  Email: tomirvine@aol.com
%
    function[f,fx,fxx,fy,fyy,fxy]= pfunction(xx,yy,a,b,r,pia,pib); 
    pia2=pia^2; 
    pib2=pib^2; 
    % 
    f(1)=1; 
    f(2)=cos(r(2)*xx); 
    f(3)=cos(r(3)*yy); 
    f(4)=cosh(r(4)*xx); 
    f(5)=cosh(r(5)*yy); 
    f(6)=cos(r(6)*xx)*cos(r(7)*yy); 
    % 
    fx(1)=0; 
    fx(2)=-pia*r(2)*sin(r(2)*xx); 
    fx(3)=0; 
    fx(4)=pia*r(4)*sinh(r(4)*xx); 
    fx(5)=0; 
    fx(6)=-pia*r(6)*sin(r(6)*xx)*cos(r(7)*yy); 
    % 
    fxx(1)=0; 
    fxx(2)=-pia2*r(2)^2*cos(r(2)*xx); 
    fxx(3)=0; 
    fxx(4)=pia2*r(4)^2*cosh(r(4)*xx); 
    fxx(5)=0; 
    fxx(6)=-pia2*r(6)^2*cos(r(6)*xx)*cos(r(7)*yy); 
    % 
    fy(1)=0; 
    fy(2)=0; 
    fy(3)=-pib*r(3)*sin(r(3)*yy); 
    fy(4)=0; 
    fy(5)=pib*r(5)*sinh(r(5)*yy); 
    fy(6)=-pib*r(7)*cos(r(6)*xx)*sin(r(7)*yy); 
    % 
    fyy(1)=0; 
    fyy(2)=0; 
    fyy(3)=-pib2*r(3)^2*cos(r(3)*yy); 
    fyy(4)=0; 
    fyy(5)=pib2*r(5)^2*cosh(r(5)*yy); 
    fyy(6)=-pib2*r(7)^2*cos(r(6)*xx)*cos(r(7)*yy); 
    % 
    fxy(1)=0; 
    fxy(2)=0; 
    fxy(3)=0; 
    fxy(4)=0; 
    fxy(5)=0; 
    fxy(6)=pia*pib*r(6)*r(7)*sin(r(6)*xx)*sin(r(7)*yy); 