Spherical harmonics. Almost everyhing is done. Missing triangle area calcul...

Samo Penic

2012-06-07 523bf18206f550a315c6c17e5a0a253381b0f8bf

Spherical harmonics. Almost everyhing is done. Missing triangle area calculation when vertex is moved or bond is flipped. Also missing volume calculation on vertex move or bondflip. Calculation of co coefficient is not done completely yet. Problems are in numbering the coefficients. Newly added data structure ts_spharm is referenced from ts_vesicle. Missing function for initialization and freeing the memory of that datastructure -- but that memory is already used by some functions.

 aclocal.m4

@@ -1,7 +1,8 @@
# generated automatically by aclocal 1.11.1 -*- Autoconf -*-
# generated automatically by aclocal 1.11.3 -*- Autoconf -*-

# Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004,
# 2005, 2006, 2007, 2008, 2009  Free Software Foundation, Inc.
# 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation,
# Inc.
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.
@@ -13,17 +14,20 @@

m4_ifndef([AC_AUTOCONF_VERSION],
  [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.67],,
[m4_warning([this file was generated for autoconf 2.67.
m4_if(m4_defn([AC_AUTOCONF_VERSION]), [2.68],,
[m4_warning([this file was generated for autoconf 2.68.
You have another version of autoconf.  It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically `autoreconf'.])])

# Copyright (C) 2002, 2003, 2005, 2006, 2007, 2008  Free Software Foundation, Inc.
# Copyright (C) 2002, 2003, 2005, 2006, 2007, 2008, 2011 Free Software
# Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 1

# AM_AUTOMAKE_VERSION(VERSION)
# ----------------------------
@@ -34,7 +38,7 @@
[am__api_version='1.11'
dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to
dnl require some minimum version.  Point them to the right macro.
m4_if([$1], [1.11.1], [],
m4_if([$1], [1.11.3], [],
      [AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])

@@ -50,18 +54,20 @@
# Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced.
# This function is AC_REQUIREd by AM_INIT_AUTOMAKE.
AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION],
[AM_AUTOMAKE_VERSION([1.11.1])dnl
[AM_AUTOMAKE_VERSION([1.11.3])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
  [m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])

# AM_AUX_DIR_EXPAND                                         -*- Autoconf -*-

# Copyright (C) 2001, 2003, 2005  Free Software Foundation, Inc.
# Copyright (C) 2001, 2003, 2005, 2011 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 1

# For projects using AC_CONFIG_AUX_DIR([foo]), Autoconf sets
# $ac_aux_dir to `$srcdir/foo'.  In other projects, it is set to
@@ -144,14 +150,14 @@
Usually this means the macro was only invoked conditionally.]])
fi])])

# Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009
# Free Software Foundation, Inc.
# Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009,
# 2010, 2011 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 10
# serial 12

# There are a few dirty hacks below to avoid letting `AC_PROG_CC' be
# written in clear, in which case automake, when reading aclocal.m4,
@@ -191,6 +197,7 @@
  # instance it was reported that on HP-UX the gcc test will end up
  # making a dummy file named `D' -- because `-MD' means `put the output
  # in D'.
  rm -rf conftest.dir
  mkdir conftest.dir
  # Copy depcomp to subdir because otherwise we won't find it if we're
  # using a relative directory.
@@ -255,7 +262,7 @@
    break
      fi
      ;;
    msvisualcpp | msvcmsys)
    msvc7 | msvc7msys | msvisualcpp | msvcmsys)
      # This compiler won't grok `-c -o', but also, the minuso test has
      # not run yet.  These depmodes are late enough in the game, and
      # so weak that their functioning should not be impacted.
@@ -320,10 +327,13 @@
if test "x$enable_dependency_tracking" != xno; then
  am_depcomp="$ac_aux_dir/depcomp"
  AMDEPBACKSLASH='\'
  am__nodep='_no'
fi
AM_CONDITIONAL([AMDEP], [test "x$enable_dependency_tracking" != xno])
AC_SUBST([AMDEPBACKSLASH])dnl
_AM_SUBST_NOTMAKE([AMDEPBACKSLASH])dnl
AC_SUBST([am__nodep])dnl
_AM_SUBST_NOTMAKE([am__nodep])dnl
])

# Generate code to set up dependency tracking.              -*- Autoconf -*-
@@ -545,11 +555,14 @@
done
echo "timestamp for $_am_arg" >`AS_DIRNAME(["$_am_arg"])`/stamp-h[]$_am_stamp_count])

# Copyright (C) 2001, 2003, 2005, 2008  Free Software Foundation, Inc.
# Copyright (C) 2001, 2003, 2005, 2008, 2011 Free Software Foundation,
# Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 1

# AM_PROG_INSTALL_SH
# ------------------
@@ -682,11 +695,14 @@
fi
])

# Copyright (C) 2003, 2004, 2005, 2006  Free Software Foundation, Inc.
# Copyright (C) 2003, 2004, 2005, 2006, 2011 Free Software Foundation,
# Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 1

# AM_PROG_MKDIR_P
# ---------------
@@ -710,13 +726,14 @@

# Helper functions for option handling.                     -*- Autoconf -*-

# Copyright (C) 2001, 2002, 2003, 2005, 2008  Free Software Foundation, Inc.
# Copyright (C) 2001, 2002, 2003, 2005, 2008, 2010 Free Software
# Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 4
# serial 5

# _AM_MANGLE_OPTION(NAME)
# -----------------------
@@ -724,13 +741,13 @@
[[_AM_OPTION_]m4_bpatsubst($1, [[^a-zA-Z0-9_]], [_])])

# _AM_SET_OPTION(NAME)
# ------------------------------
# --------------------
# Set option NAME.  Presently that only means defining a flag for this option.
AC_DEFUN([_AM_SET_OPTION],
[m4_define(_AM_MANGLE_OPTION([$1]), 1)])

# _AM_SET_OPTIONS(OPTIONS)
# ----------------------------------
# ------------------------
# OPTIONS is a space-separated list of Automake options.
AC_DEFUN([_AM_SET_OPTIONS],
[m4_foreach_w([_AM_Option], [$1], [_AM_SET_OPTION(_AM_Option)])])
@@ -806,11 +823,13 @@
fi
AC_MSG_RESULT(yes)])

# Copyright (C) 2001, 2003, 2005  Free Software Foundation, Inc.
# Copyright (C) 2001, 2003, 2005, 2011 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 1

# AM_PROG_INSTALL_STRIP
# ---------------------
@@ -834,13 +853,13 @@
INSTALL_STRIP_PROGRAM="\$(install_sh) -c -s"
AC_SUBST([INSTALL_STRIP_PROGRAM])])

# Copyright (C) 2006, 2008  Free Software Foundation, Inc.
# Copyright (C) 2006, 2008, 2010 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
# with or without modifications, as long as this notice is preserved.

# serial 2
# serial 3

# _AM_SUBST_NOTMAKE(VARIABLE)
# ---------------------------
@@ -849,13 +868,13 @@
AC_DEFUN([_AM_SUBST_NOTMAKE])

# AM_SUBST_NOTMAKE(VARIABLE)
# ---------------------------
# --------------------------
# Public sister of _AM_SUBST_NOTMAKE.
AC_DEFUN([AM_SUBST_NOTMAKE], [_AM_SUBST_NOTMAKE($@)])

# Check how to create a tarball.                            -*- Autoconf -*-

# Copyright (C) 2004, 2005  Free Software Foundation, Inc.
# Copyright (C) 2004, 2005, 2012 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
@@ -877,10 +896,11 @@
# a tarball read from stdin.
#     $(am__untar) < result.tar
AC_DEFUN([_AM_PROG_TAR],
[# Always define AMTAR for backward compatibility.
AM_MISSING_PROG([AMTAR], [tar])
[# Always define AMTAR for backward compatibility.  Yes, it's still used
# in the wild :-(  We should find a proper way to deprecate it ...
AC_SUBST([AMTAR], ['$${TAR-tar}'])
m4_if([$1], [v7],
     [am__tar='${AMTAR} chof - "$$tardir"'; am__untar='${AMTAR} xf -'],
     [am__tar='$${TAR-tar} chof - "$$tardir"' am__untar='$${TAR-tar} xf -'],
     [m4_case([$1], [ustar],, [pax],,
              [m4_fatal([Unknown tar format])])
AC_MSG_CHECKING([how to create a $1 tar archive])

 src/general.h

@@ -110,6 +110,23 @@

/* STRUCTURES */


/** @brief Data structure for keeping the coordinates in selected coordinate
 * system
 */
#define TS_COORD_CARTESIAN 0
#define TS_COORD_SPHERICAL 1
#define TS_COORD_CYLINDRICAL 2

typedef struct {
    ts_double e1;
    ts_double e2;
    ts_double e3;
    ts_uint coord_type;
} ts_coord;



/** @brief Data structure of all data connected to a vertex
 *
 *  ts_vertex holds the data for one single point (bead, vertex). To understand how to use it
@@ -134,6 +151,9 @@
        ts_double xk;
        ts_double c;
        ts_uint id;
        ts_double projArea;
        ts_double relR;
        ts_double solAngle;
};
typedef struct ts_vertex ts_vertex;

@@ -166,6 +186,7 @@
    ts_double xnorm;
    ts_double ynorm;
    ts_double znorm;
    ts_double area; // firstly needed for sh.c
};
typedef struct ts_triangle ts_triangle;

@@ -193,6 +214,16 @@


typedef struct {
    ts_uint l;
    ts_uint i;
    ts_double ***Ylmi;
    ts_double **ulm;
    ts_uint **co;
} ts_spharm;



typedef struct {
    ts_vertex_list *vlist;
    ts_bond_list *blist;
    ts_triangle_list *tlist;
@@ -202,9 +233,13 @@
    ts_double dmax;
    ts_double stepsize;
    ts_double cm[3];
    ts_double volume;
    ts_spharm *sphHarmonics;
} ts_vesicle;




/* GLOBAL VARIABLES */

int quiet;

 src/sh.c

@@ -55,6 +55,20 @@
}



ts_bool precomputeShCoeff(ts_spharm *sph){
    ts_uint i,j;
    for(i=0;i<sph->l;i++){
        sph->co[i][i]=sqrt((2.0*i+1.0)/2.0/M_PI);
        for(j=0;j<i-1;j++){
            
        }
    }
    return TS_SUCCESS;

}


/*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;
@@ -88,3 +102,148 @@
    
    return K*sqrt((2.0*l+1.0)/2.0*fac1/fac2)*plgndr(l,abs(m),cos(theta));    
}


/* Function transforms coordinates from cartesian to spherical coordinates
 * (r,phi, theta). */
ts_bool *cart2sph(ts_coord *coord, ts_double x, ts_double y, ts_double z){
    coord->coord_type=TS_COORD_SPHERICAL;
#ifdef TS_DOUBLE_DOUBLE
    coord->e1=sqrt(x*x+y*y+z*z);
    if(z==0) coord->e3=M_PI/2.0;
    else coord->e3=atan(sqrt(x*x+y*y)/z);
    coord->e2=atan2(y,x);
#endif
#ifdef TS_DOUBLE_FLOAT
    coord->e1=sqrtf(x*x+y*y+z*z);
    if(z==0) coord->e3=M_PI/2.0;
    else coord->e3=atanf(sqrtf(x*x+y*y)/z);
    coord->e2=atan2f(y,x);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
    coord->e1=sqrtl(x*x+y*y+z*z);
    if(z==0) coord->e3=M_PI/2.0;
    else coord->e3=atanl(sqrtl(x*x+y*y)/z);
    coord->e2=atan2l(y,x);
#endif

    return TS_SUCCESS;
}

/* Function returns radius of the sphere with the same volume as vesicle (r0) */
ts_double getR0(ts_vesicle *vesicle){
    ts_double r0;
 #ifdef TS_DOUBLE_DOUBLE
   r0=pow(vesicle->volume*3.0/4.0/M_PI,1.0/3.0);
#endif
#ifdef TS_DOUBLE_FLOAT
   r0=powf(vesicle->volume*3.0/4.0/M_PI,1.0/3.0);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
   r0=powl(vesicle->volume*3.0/4.0/M_PI,1.0/3.0);
#endif
    return r0;
}


ts_bool preparationSh(ts_vesicle *vesicle, ts_double r0){
//TODO: before calling or during the call calculate area of each triangle! Can
//be also done after vertexmove and bondflip //
    ts_uint i,j;
    ts_vertex **vtx=vesicle->vlist->vtx;
    ts_vertex *cvtx;
    ts_triangle *ctri;
    ts_double centroid[3];
    ts_double r;
    for (i=0;  i<vesicle->vlist->n; i++){
        cvtx=vtx[i];
        //cvtx->projArea=4.0*M_PI/1447.0*(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z)/r0/r0;
        cvtx->projArea=0.0;

        /* go over all triangles that have a common vertex i */
        for(j=0; j<cvtx->tristar_no; j++){
            ctri=cvtx->tristar[j];
            centroid[0]=(ctri->vertex[0]->x + ctri->vertex[1]->x + ctri->vertex[2]->x)/3.0;
            centroid[1]=(ctri->vertex[0]->y + ctri->vertex[1]->y + ctri->vertex[2]->y)/3.0;
            centroid[2]=(ctri->vertex[0]->z + ctri->vertex[1]->z + ctri->vertex[2]->z)/3.0;
        /* calculating projArea+= area(triangle)*cos(theta) */
#ifdef TS_DOUBLE_DOUBLE
            cvtx->projArea = cvtx->projArea + ctri->area*(-centroid[0]*ctri->xnorm - centroid[1]*ctri->ynorm - centroid[2]*ctri->znorm)/ sqrt(centroid[0]*centroid[0]+centroid[1]*centroid[1]+centroid[2]*centroid[2]);
#endif
#ifdef TS_DOUBLE_FLOAT
            cvtx->projArea = cvtx->projArea + ctri->area*(-centroid[0]*ctri->xnorm - centroid[1]*ctri->ynorm - centroid[2]*ctri->znorm)/ sqrtf(centroid[0]*centroid[0]+centroid[1]*centroid[1]+centroid[2]*centroid[2]);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
            cvtx->projArea = cvtx->projArea + ctri->area*(-centroid[0]*ctri->xnorm - centroid[1]*ctri->ynorm - centroid[2]*ctri->znorm)/ sqrtl(centroid[0]*centroid[0]+centroid[1]*centroid[1]+centroid[2]*centroid[2]);
#endif
        }

    cvtx->projArea=cvtx->projArea/3.0;
        //we dont store spherical coordinates of vertex, so we have to calculate
        //r(i) at this point.
#ifdef TS_DOUBLE_DOUBLE
    r=sqrt(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
#endif
#ifdef TS_DOUBLE_FLOAT
    r=sqrtf(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
#endif
#ifdef TS_DOUBLE_LONGDOUBLE
    r=sqrtl(cvtx->x*cvtx->x+cvtx->y*cvtx->y+cvtx->z*cvtx->z);
#endif
    cvtx->relR=(r-r0)/r0;
    cvtx->solAngle=cvtx->projArea/cvtx->relR * cvtx->projArea/cvtx->relR;
    }
    return TS_SUCCESS;
}



ts_bool calculateYlmi(ts_vesicle *vesicle){
    ts_uint i,j,k;
    ts_spharm *sph=vesicle->sphHarmonics;
    ts_coord *coord=(ts_coord *)malloc(sizeof(ts_coord));
    ts_double fi, theta;
    for(k=0;k<vesicle->vlist->n;k++){
        sph->Ylmi[0][0][k]=sqrt(1.0/4.0/M_PI);
        cart2sph(coord,vesicle->vlist->vtx[k]->x, vesicle->vlist->vtx[k]->y, vesicle->vlist->vtx[k]->z);
        fi=coord->e2;
        theta=coord->e3; 
        for(i=0; i<sph->l; i++){
            for(j=0;j<i;j++){
                sph->Ylmi[i][j][k]=sph->co[i][j]*cos((j-i-1)*fi)*pow(-1,j-i-1)*plgndr(i,abs(j-i-1),cos(theta));
            }
                sph->Ylmi[i][j+1][k]=sph->co[i][j+1]*plgndr(i,0,cos(theta));
            for(j=sph->l;j<2*i;j++){
                sph->Ylmi[i][j][k]=sph->co[i][j]*sin((j-i-1)*fi)*plgndr(i,j-i-1,cos(theta));
            }
        }

    }
    free(coord);
    return TS_SUCCESS;
}



ts_bool calculateUlm(ts_vesicle *vesicle){
    ts_uint i,j,k;
    ts_vertex *cvtx;
    for(i=0;i<vesicle->sphHarmonics->l;i++){
        for(j=0;j<2*i;j++) vesicle->sphHarmonics->ulm[i][j]=0.0;
    }

//TODO: call calculateYlmi !!!


    for(k=0;k<vesicle->vlist->n; k++){
        cvtx=vesicle->vlist->vtx[k];
        for(i=0;i<vesicle->sphHarmonics->l;i++){
            for(j=0;j<2*i;j++){
                vesicle->sphHarmonics->ulm[i][j]+= cvtx->solAngle*cvtx->relR*vesicle->sphHarmonics->Ylmi[i][j][k];
            }

        }
    }

    return TS_SUCCESS;
}

 src/vesicle.c

@@ -36,3 +36,21 @@
    free(vesicle);
    return TS_SUCCESS;
}

ts_bool vesicle_volume(ts_vesicle *vesicle){
    ts_double volume;
    ts_double vol;
    ts_uint i;
    ts_triangle **tria=vesicle->tlist->tria;
    volume=0;
    for(i=0; i<vesicle->tlist->n;i++){
        vol=(tria[i]->vertex[0]->x+ tria[i]->vertex[1]->x + tria[i]->vertex[2]->x) * tria[i]->xnorm + 
       (tria[i]->vertex[0]->y+ tria[i]->vertex[1]->y + tria[i]->vertex[2]->y) * tria[i]->ynorm + 
    (tria[i]->vertex[0]->z+ tria[i]->vertex[1]->z + tria[i]->vertex[2]->z) *
tria[i]->znorm;
    volume=volume-vol/18.0;
    }

    vesicle->volume=volume;
    return TS_SUCCESS;
}