c#include <fintrf.h>
C     gateway subroutine

      subroutine mexfunction(nlhs, plhs, nrhs, prhs)
      implicit none
C-----------------------------------------------------------------------
C     (integer) Replace integer by integer*8 on 64-bit platforms
C
      integer plhs(*), prhs(*)
      integer mxGetPr, mxCreateDoubleMatrix, mxIsChar, mxGetString
      integer mxGetM, mxGetN

      integer i_r, i_theta, i_phi, i_xmaxdeg
      integer i_A_r, i_A_theta, i_A_phi
C-----------------------------------------------------------------------
C
      integer nlhs, nrhs
      integer m_in, n_in, size2, size3
      integer N_data,maxdeg,nkoeff
      real*8 xmaxdeg
      character*3 source
      integer  status, strlen

      if(nrhs .ne. 5) then
         call mexErrMsgTxt('5 inputs required.')
      endif
      if(nlhs .ne. 3) then
         call mexErrMsgTxt('Three outputs required.')
      endif
         
      m_in = mxGetM(prhs(1))
      n_in = mxGetN(prhs(1))
      N_data = m_in * n_in
      i_r = mxGetPr(prhs(1))

      m_in = mxGetM(prhs(2))
      n_in = mxGetN(prhs(2))
      size2 = m_in * n_in
      i_theta = mxGetPr(prhs(2))

      m_in = mxGetM(prhs(3))
      n_in = mxGetN(prhs(3))
      size3 = m_in * n_in
      i_phi = mxGetPr(prhs(3))

      if((N_data.ne.size2).or.(N_data.ne.size3)) print*,
     *            'r, theta, phi must have the same length'

      m_in = mxGetM(prhs(4))
      n_in = mxGetN(prhs(4))
      if(m_in.ne.1) print*,' 4th argument not a scalar'
      if(n_in.ne.1) print*,' 4th argument not a scalar'
      i_xmaxdeg = mxGetPr(prhs(4))
      call mxCopyPtrToReal8(i_xmaxdeg, xmaxdeg, 1)
      maxdeg = xmaxdeg

C     The input  must be a string.
      if(mxIsChar(prhs(5)) .ne. 1) then
         call mexErrMsgTxt('Input must be a string.')
C     The input must be a row vector.
      elseif (mxGetM(prhs(5)) .ne. 1) then
         call mexErrMsgTxt('Input must be a row vector.')
      endif
    
C     Get the length of the input string.
      strlen = mxGetM(prhs(5))*mxGetN(prhs(5))
C     Get the string contents (dereference the input integer).
      status = mxGetString(prhs(5), source, 3)
C     Check if mxGetString is successful.
      if (status .ne. 0) then 
         call mexErrMsgTxt('Filename length must be less than 3.')
      endif

      nkoeff=maxdeg*(maxdeg+2)
      plhs(1) = mxCreateDoubleMatrix(N_data, nkoeff, 0)
      i_A_r = mxGetPr(plhs(1))
      plhs(2) = mxCreateDoubleMatrix(N_data, nkoeff, 0)
      i_A_theta = mxGetPr(plhs(2))
      plhs(3) = mxCreateDoubleMatrix(N_data, nkoeff, 0)
      i_A_phi = mxGetPr(plhs(3))

C     call the computational routine
      call make_dm(%val(i_A_r), %val(i_A_theta), %val(i_A_phi), 
     *      %val(i_r), %val(i_theta), %val(i_phi), 
     *      N_data, %val(i_xmaxdeg), source)

      return
      end
      
C     computational subroutine
      subroutine make_dm(A_r, A_theta, A_phi, 
     *     r, theta, phi, 
     *     N_data, xmaxdeg, source)
	implicit real*8 (a-h,o-z)
      real*8 r(*), theta(*), phi(*)
      real*8 A_r(N_data,*), A_theta(N_data,*), A_phi(N_data,*)

	CHARACTER*3 source
	integer kn,maxdeg
	real*8 a,pi,rad
	parameter(kn=750,nin=kn*(kn+2))
	parameter (pi=3.14159265358979d0,rad=pi/180.d0)

c other declarations
	integer nu(nin),mu(nin),l,n,m,lmax
	real*8 pd(0:kn,0:kn+1)
	real*8 a__r(kn),cml(0:kn),sml(0:kn)
	real*8 sinb,cosb

	cml(0)=1.		! Initialisation cos(0.)
	sml(0)=0.		! Initialisation sin(0.)

      maxdeg = xmaxdeg
      if(maxdeg .gt. kn) call mexErrMsgTxt('Maxdeg > kn in design_SHA')

	l=1			! List of coefficients to be determined
	do  m=0,maxdeg
	   do  n=max(1,m),maxdeg
	      nu(l)=n
	      mu(l)=m
	      l=l+1
	   enddo
	enddo
	lmax=l-1

      do i=1, N_data
	if(source(1:1).eq.'e'.or.source(1:1).eq.'E') then ! external sources
	do  n=1,maxdeg
	   cml(n)=cos(n*phi(i))
	   sml(n)=sin(n*phi(i))
	   a__r(n)=r(i)**(n-1)
	enddo
	elseif(source(1:1).eq.'t'.or.source(1:1).eq.'T') then ! toroidal sources
	   cml(n)=cos(n*phi(i))
	   sml(n)=sin(n*phi(i))
	   a__r(n)=r(i)**(-1)
      else ! internal sources as default
	do  n=1,maxdeg
	   cml(n)=cos(n*phi(i))
	   sml(n)=sin(n*phi(i))
	   a__r(n)=r(i)**(-(n+2))
	enddo
	endif
	cosb=cos(theta(i))
	sinb=sin(theta(i))
      call pnmas (cosb,maxdeg,kn,pd)
c
c construct A_r, A_theta, A_phi
c
	l=0
	if(source(1:1).eq.'e'.or.source(1:1).eq.'E') then
	do n=1,maxdeg
	   do m=0,n
	      l=l+1
	      A_r(i,l)=-n*pd(n,m)*cml(m)*a__r(n)
	      A_theta(i,l)=-pd(m,n+1)*cml(m)*a__r(n)
	      A_phi(i,l)=pd(n,m)*sml(m)*m/sinb*a__r(n)
	      if(m.gt.0) then
		 l=l+1
		 A_r(i,l)=-n*pd(n,m)*sml(m)*a__r(n)
		 A_theta(i,l)=-pd(m,n+1)*sml(m)*a__r(n)
		 A_phi(i,l)=-pd(n,m)*cml(m)*a__r(n)*m/sinb
	      endif
	   enddo
	enddo
	elseif(source(1:1).eq.'t'.or.source(1:1).eq.'T') then
	do n=1,maxdeg
	   do m=0,n
	      l=l+1
	      A_r(i,l)=0
	      A_theta(i,l)=pd(n,m)*sml(m)*m/sinb*a__r(n)
	      A_phi(i,l)  =-pd(m,n+1)*cml(m)*a__r(n)
	      if(m.gt.0) then
		 l=l+1
		 A_r(i,l)=0
		 A_theta(i,l)=-pd(n,m)*cml(m)*a__r(n)*m/sinb
		 A_phi(i,l)  =-pd(m,n+1)*sml(m)*a__r(n)
	      endif
	   enddo
	enddo
	else ! internal sources as default
	do n=1,maxdeg
	   do m=0,n
	      l=l+1
	      A_r(i,l)=+(n+1.)*pd(n,m)*cml(m)*a__r(n)
	      A_theta(i,l)=-pd(m,n+1)*cml(m)*a__r(n)
	      A_phi(i,l)=pd(n,m)*sml(m)*m/sinb*a__r(n)
	      if(m.gt.0) then
		 l=l+1
		 A_r(i,l)=+(n+1.)*pd(n,m)*sml(m)*a__r(n)
		 A_theta(i,l)=-pd(m,n+1)*sml(m)*a__r(n)
		 A_phi(i,l)=-pd(n,m)*cml(m)*a__r(n)*m/sinb
	      endif
	   enddo
	enddo
	endif
      enddo ! loop over N_data

      return
      end

      SUBROUTINE PNMAS(X,Maxdeg,Kn,Pd)
	implicit real*8 (a-h,o-z)
C
C   Calculation of Legendre function P(N,M) (Schmidt normalized)
C   and its derivative dP(N,M) vrt. theta.
C
C   Input:  X     real*8      cos(theta)
C           MAXDEG  integer*4   maximum degree
C           KN    integer*4   dimension of array containing Legendre functions
C                             is PD(0:KN, 0:KN+1)
C   Output: PD    real*8      array PD(0:KN, 0:KN+1) with Legendre functions
C
C   P(n,m) ==> PD(n,m) and dP(n,m) ==> PD(m,n+1)
C   P(N,M) are the associated Legendre functions (Schmidt-normalized)
C
C   March 1997, Nils Olsen, NBI Copenhagen.
C
      INTEGER Kn , Maxdeg , n , m , l
      REAL*8 Pd(0:Kn,0:Kn+1) , X , sinth , pmm , e , d
C
      sinth = SQRT((1.D0-X)*(1.D0+X))
      Pd(0,0) = 1.D0 ! starting values
      Pd(1,1) = sinth
c Recursion relations after Langel "The Main Field" (1987),
c eq. (27) and Table 2 (p. 256)
      DO 100 m = 0 , Maxdeg - 1
         pmm = SQRT(m+m+1.D0)*Pd(m,m)
         Pd(m+1,m) = X*pmm
         IF ( m.GT.0 ) Pd(m+1,m+1) = sinth*pmm/SQRT(m+m+2.D0)
         DO 50 n = m + 2 , Maxdeg
            d = n*n - m*m
            e = n + n - 1.D0
            Pd(n,m) = (e*X*Pd(n-1,m)-SQRT(d-e)*Pd(n-2,m))/SQRT(d)
 50      CONTINUE
 100  CONTINUE
C
C     dP(N,M)=PD(M,N+1) is the derivative of P(N,M) vrt. theta
C
      Pd(0,2) = -Pd(1,1)
      Pd(1,2) =  Pd(1,0)
      DO 200 n = 2 , Maxdeg
         l = n + 1
         Pd(0,l) = -SQRT(.5D0*(n*n+n))*Pd(n,1)
         Pd(1,l) = .5D0*(SQRT(2.D0*(n*n+n))*Pd(n,0)-SQRT((n*n+n-2.D0))
     &             *Pd(n,2))
         DO 150 m = 2 , n - 1
            Pd(m,l) = .5D0*(SQRT((n+m)*(n-m+1.D0))*Pd(n,m-1)
     &                -SQRT((n+m+1.D0)*(n-m))*Pd(n,m+1))
 150     CONTINUE
         Pd(n,l) = .5D0*SQRT(2.D0*n)*Pd(n,n-1)
 200  CONTINUE
      END


