proj/globalhs/libs/projection/setplm2.c

// $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/projection/setplm2.c,v 1.1 2013/04/28 07:46:05 tplarson Exp $

// this function is a port of the one described below,
// and adds calcuation of derivative

/*
c
c  plm1.f                            ~soi/(version)/src/functions/plm1.f
c
c  Description:
c     source: Jesper Schou
c     sets plm(i,indx(l))=plm(x(i),l) for l=lmin,lmax and i=1,nx
c     the plm's are normalized such that int from -1 to 1 plm sup 2 dx=1
c     all 0<=m<=l<=1000 tested with nx=5001, max err always at m=0
c     all 0<=m<=10 and m<=l<=2500 tested with nx=12501
c     every 50 l and m with 0<=m<=l<=3000 tested with nx=10001
c     fails at l \approx 1800 and m=700 because of underflow in setting
c     up P_m^m
c     if lmin>m then the Plm's for l=lmin-2 (if applicable) and
c     l=lmin-1 must be set on input
c
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define minval(x,y)     (((x) < (y)) ? (x) : (y))
#define maxval(x,y)     (((x) < (y)) ? (y) : (x))

#define MAXIMUML 1800

char *cvsinfo_setplm2 = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/libs/projection/setplm2.c,v 1.1 2013/04/28 07:46:05 tplarson Exp $";

setplm2(int lmin,int lmax,int m,long nx,int *indx,double *x,long nplm,double *plm,double *dplm)
{
  double eps=1.0e-12;
  double *x1=malloc(nx*sizeof(double));
  double *x2=malloc(nx*sizeof(double));
  double *x3=malloc(nx*sizeof(double));
  if (lmax > MAXIMUML)
  {
    fprintf(stderr,"setplm2 called with lmax larger than %d\n", MAXIMUML);
    return 1;
  }

  if (lmin < m)
    lmin=m;

  long i,j;
  for (i=0;i<nx;i++)
  {
    x1[i]=minval(1.0-eps,maxval(-1.0+eps,x[i]));
    x2[i]=1.0/(x1[i]*x1[i]-1.0);
    x3[i]=x1[i]*x2[i];
  }

  int im,im1;
  double c;
  if (lmin <= m+1)
  {
    im=indx[m];
    c=sqrt((2*m+1)/2.0);
    for (i=1;i<=m;i++)
      c*=-sqrt(1-0.5/i);
  }

  if (lmin == m)
  {
    for (i=0;i<nx;i++)
    {
//      plm[im][i]=c;
      plm[im*nplm + i]=c;
      for (j=1;j<=m;j++)
//        plm[im][i]*=sqrt(1.0 - x1[i]*x1[i]);
        plm[im*nplm + i]*=sqrt(1.0 - x1[i]*x1[i]);
    }
  }

  if (lmax > m && lmin <= m+1)
  {
    im1=indx[m+1];
    c=sqrt(2*m + 3.0);
    for (i=0;i<nx;i++)
//      plm[im1][i]=x1[i]*c*plm[im][i];
      plm[im1*nplm + i]=x1[i]*c*plm[im*nplm + i];
  }

  long m2,l,l2,il0,il1,il2;
  double c1,c2;
  m2=m*m;
  for (l=maxval(lmin,m+2);l<=lmax;l++)
  {
    l2=l*l;
    il0=indx[l];
    il1=indx[l-1];
    il2=indx[l-2];
    c1=sqrt((4*l2-1.0)/(l2-m2));
    c2=sqrt(((2*l+1.0)*(l+m-1)*(l-m-1))/((2*l-3.0)*(l2-m2)));
    for (i=0;i<nx;i++)
//      plm[il0][i]=c1*x1[i]*plm[il1][i]-c2*plm[il2][i];
      plm[il0*nplm + i]=c1*x1[i]*plm[il1*nplm + i]-c2*plm[il2*nplm + i];
  }

  if (dplm != NULL)
  {
    im=indx[m];
    for (i=0;i<nx;i++)
//      dplm[im][i]=m*x3[i]*plm[im][i];
      dplm[im*nplm + i]=m*x3[i]*plm[im*nplm + i];
    for (l=m+1;l<=lmax;l++)
    {
      il0=indx[l];
      il1=indx[l-1];
      c=sqrt((2.0*l+1.0)*(l*l-m2)/(2.0*l-1));
      for (i=0;i<nx;i++)
//        dplm[il0][i]=l*x3[i]*plm[il0][i]-c*x2[i]*plm[il1][i];
        dplm[il0*nplm + i]=l*x3[i]*plm[il0*nplm + i]-c*x2[i]*plm[il1*nplm + i];

    }
  }

  free(x1);
  free(x2);
  free(x3);

}

---