Added spherical harmonics routine.I've got a feeling that some harmonics ar...

Samo Penic

2011-03-19 88f45143117f3ae973be6e3ea5dac8063a7b8ece

Added spherical harmonics routine.I've got a feeling that some harmonics are not calculated correctly!

 src/Makefile.am

@@ -1,5 +1,8 @@
trisurfdir=../
trisurf_PROGRAMS = trisurf
trisurf_SOURCES = general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c frame.c energy.c timestep.c vertexmove.c main.c
trisurf_LDFLAGS = -lm -lconfuse
trisurf_CFLAGS = -Wall -Werror
#trisurf_LDFLAGS = -lm -lconfuse
shdiscoverdir=../
shdiscover_PROGRAMS= shdiscover
shdiscover_SOURCES= general.c vertex.c bond.c triangle.c cell.c vesicle.c initial_distribution.c io.c energy.c sh.c shdiscover.c
AM_CFLAGS = -Wall -Werror

 src/sh.c

New file
@@ -0,0 +1,90 @@
#include<math.h>
#include<stdlib.h>
#include "general.h"
#include "sh.h"

/* Gives you legendre polynomials. Taken from NR, p. 254 */
ts_double plgndr(ts_int l, ts_int m, ts_float x){
    ts_double fact, pll, pmm, pmmp1, somx2;
    ts_int i,ll;

#ifdef TS_DOUBLE_DOUBLE
    if(m<0 || m>l || fabs(x)>1.0)
        fatal("Bad arguments in routine plgndr",1);
#endif
#ifdef TS_DOUBLE_FLOAT
    if(m<0 || m>l || fabsf(x)>1.0)
        fatal("Bad arguments in routine plgndr",1);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
    if(m<0 || m>l || fabsl(x)>1.0)
        fatal("Bad arguments in routine plgndr",1);
#endif
    pmm=1.0;
    if (m>0) {
#ifdef TS_DOUBLE_DOUBLE
        somx2=sqrt((1.0-x)*(1.0+x));
#endif
#ifdef TS_DOUBLE_FLOAT
        somx2=sqrtf((1.0-x)*(1.0+x));
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
        somx2=sqrtl((1.0-x)*(1.0+x));
#endif
        fact=1.0;
        for (i=1; i<=m;i++){
            pmm *= -fact*somx2;
            fact +=2.0;
        }
    }

    if (l == m) return pmm;
    else {
        pmmp1=x*(2*m+1)*pmm;
        if(l==(m+1)) return(pmmp1);
        else {
            pll=0; /* so it can not be uninitialized */
            for(ll=m+2;ll<=l;ll++){
                pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
                pmm=pmmp1;
                pmmp1=pll;
            }
            return(pll);
        }
    }
}


/*Computes Y(l,m,theta,fi) (Miha's definition that is different from common definition for  factor srqt(1/(2*pi)) */
ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi){
    ts_double fac1, fac2, K;
    int i;

    if(l<0 || m>l || m<-l)
        fatal("Error using shY function!",1);

    fac1=1.0;
    for(i=1; i<=l-m;i++){
        fac1 *= i;
    }
    fac2=1.0;
    for(i=1; i<=l+m;i++){
        fac2 *= i;
    }

    if(m==0){
        K=sqrt(1.0/(2.0*M_PI));
    }
    else if (m>0) {
        K=sqrt(1.0/(M_PI))*cos(m*fi);
    } 
    else {
        //K=pow(-1.0,abs(m))*sqrt(1.0/(2.0*M_PI))*cos(m*fi);
        if(abs(m)%2==0)
        K=sqrt(1.0/(M_PI))*sin(m*fi);
        else
        K=-sqrt(1.0/(M_PI))*sin(m*fi);
    }
	
    return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));	
}

 src/sh.h

New file
@@ -0,0 +1,9 @@
#ifndef _H_SH
#define _H_SH
#include "general.h"

ts_double plgndr(ts_int l, ts_int m, ts_float x);
ts_double shY(ts_int l,ts_int m,ts_double theta,ts_double fi);


#endif

 src/shdiscover.c

New file
@@ -0,0 +1,80 @@
#include "general.h"
#include "vertex.h"
#include "initial_distribution.h"
#include "io.h"
#include "vesicle.h"
#include "sh.h"
#include "frame.c"
#include <math.h>
#include <stdlib.h>
int main(int argc, char *argv[]){

ts_fprintf(stdout,"SHdiscover was called with %d coefficients!\n",argc-1);
ts_uint n,i,j,l;
ts_int m;
ts_double fi,theta,r,Y;
ts_vesicle *vesicle=initial_distribution_dipyramid(17,60,60,60,0.15);
ts_vertex_list *vlist=vesicle->vlist;
centermass(vesicle);
ts_fprintf(stdout,"Vesicle has a CenterMass in %f,%f,%f\n",vesicle->cm[0],vesicle->cm[1], vesicle->cm[2]);

n=vlist->n;

ts_fprintf(stdout,"Tests\n");
ts_fprintf(stdout,"P(0,0,0.5)=%f (%f)\n",plgndr(0,0,0.5),1.0);
ts_fprintf(stdout,"P(1,0,0.5)=%f (%f)\n",plgndr(1,0,0.5),0.5);
ts_fprintf(stdout,"P(2,0,0.5)=%f (%f)\n",plgndr(2,0,0.5),0.5*(3*0.5*0.5-1));

ts_fprintf(stdout,"Y(0,0,pi/6,pi/4)=%f (%f)\n",shY(0,0,M_PI/6,M_PI/4),sqrt(1/(4*M_PI)));
ts_fprintf(stdout,"Y(1,0,pi/6,pi/4)=%f (%f)\n",shY(1,0,M_PI/6,M_PI/4),sqrt(3/(4*M_PI))*cos(M_PI/6));
ts_fprintf(stdout,"Y(1,0,4*pi/6,6*pi/4)=%f (%f)\n",shY(1,0,4*M_PI/6,6*M_PI/4),sqrt(3/(4*M_PI))*cos(4*M_PI/6));
ts_fprintf(stdout,"Y(1,1,pi/6,pi/4)=%f (%f)\n",shY(1,1,M_PI/6,M_PI/4),-sqrt(3/(8*M_PI))*sin(M_PI/6)*sin(M_PI/4));
ts_fprintf(stdout,"Y(2,0,pi/6,pi/4)=%f (%f)\n",shY(2,0,M_PI/6,M_PI/4),sqrt(5/(4*M_PI))*(3.0/2.0*cos(M_PI/6)*cos(M_PI/6)-1.0/2.0));
ts_fprintf(stdout,"Y(2,-2,pi/6,pi/4)=%f (0)\n",shY(2,-2,M_PI/6,M_PI/4));
ts_fprintf(stdout,"Y(2,2,pi/6,pi/4)=%f (0)\n",shY(2,2,M_PI/6,M_PI/4));
	
    for(j=1;j<argc;j++){
        l=(int)sqrt(j-1); /* determine l from dataline */
        m=j-1-l*(l+1); /* determine m from dataline */
        ts_fprintf(stdout,"l=%d, m=%d, u=%s\n",l,m,argv[j]);
    }

/*we calculate new position of each vertex of vesicle */
for(i=0;i<n;i++){
    fi=atan2(vlist->vtx[i]->data->y, vlist->vtx[i]->data->x);
/*    theta=atan2(
        sqrt(vlist->vtx[i]->data->x*vlist->vtx[i]->data->x + 
        vlist->vtx[i]->data->y*vlist->vtx[i]->data->y),
        vlist->vtx[i]->data->z 
        ); */
    theta=acos(
        vlist->vtx[i]->data->z /
        sqrt(vlist->vtx[i]->data->x*vlist->vtx[i]->data->x + 
        vlist->vtx[i]->data->y*vlist->vtx[i]->data->y+
        vlist->vtx[i]->data->z*vlist->vtx[i]->data->z)

        );



    r=0.0;
    for(j=1;j<argc;j++){
        l=(int)sqrt(j-1); /* determine l from dataline */
        m=j-1-l*(l+1); /* determine m from dataline */
        Y=shY(l,m,theta,fi);
        r+=fabs(atof(argv[j])*Y);
        /*ts_fprintf(stdout,"l=%d, m=%d, u=%s\n",l,m,argv[j]);*/
    }

    vlist->vtx[i]->data->z=fabs(r)*cos(theta);
    vlist->vtx[i]->data->x=fabs(r)*sin(theta)*cos(fi);
    vlist->vtx[i]->data->y=fabs(r)*sin(theta)*sin(fi);
}

write_vertex_xml_file(vesicle,0);
write_master_xml_file("test.pvd");


vesicle_free(vesicle);
return 0;
}