solidangle.c

Go to the documentation of this file.

/*
    This file is part of magpar.

    Copyright (C) 2002-2010 Werner Scholz

    www:   http://www.magpar.net/
    email: magpar(at)magpar.net

    magpar is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    magpar is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with magpar; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

static char const Id[] = "$Id: solidangle.c 2962 2010-02-04 19:50:44Z scholz $\n\n";
static char const Td[] = "$Today: " __FILE__ "  " __DATE__ "  "  __TIME__  " $\n\n";

#include "util.h"

int PointFromPlane(PetscReal *x, PetscReal *v1, PetscReal *v2, PetscReal *v3, PetscReal *d)
{
  PetscReal  ab[ND],ac[ND];
  PetscReal  n[ND];

  /* calculate edge vectors */
  my_dcopy(ND,v1,1,ab,1);
  my_daxpy(ND,-1.0,v2,1,ab,1);
  my_dcopy(ND,v1,1,ac,1);
  my_daxpy(ND,-1.0,v3,1,ac,1);

  /* calculate normal vector */
  douter(ND,ab,ac,n);

  /* calculate distance */
  *d=my_ddot(ND,x,1,n,1)-my_ddot(ND,v1,1,n,1);

  return(0);
}

int SolidAngle(PetscReal *x, PetscReal *v1, PetscReal *v2, PetscReal *v3, PetscReal *omega)
{
  /* http://en.wikipedia.org/wiki/Solid_angle */

  PetscReal d;
  PointFromPlane(x,v1,v2,v3,&d);
  if (PetscAbsReal(d)<D_EPS) {
    *omega=0.0;
    return(0);
  }

  PetscReal  t_ea[ND],t_eb[ND],t_ec[ND];
  PetscReal  t_nab[ND],t_nbc[ND],t_nca[ND];
  PetscReal  t_norm;

  /* calculate edge vectors */
  my_dcopy(ND,v1,1,t_ea,1);
  my_daxpy(ND,-1.0,x,1,t_ea,1);
  my_dcopy(ND,v2,1,t_eb,1);
  my_daxpy(ND,-1.0,x,1,t_eb,1);
  my_dcopy(ND,v3,1,t_ec,1);
  my_daxpy(ND,-1.0,x,1,t_ec,1);

  /* calculate normal vectors */
  douter(ND,t_ea,t_eb,t_nab);
  douter(ND,t_eb,t_ec,t_nbc);
  douter(ND,t_ec,t_ea,t_nca);

  /* normalize vectors */

  t_norm=my_dnrm2(ND,t_nab,1);
  if (t_norm < D_EPS) {
    *omega=0.0;
    return(0);
  }
  my_dscal(ND,1.0/t_norm,t_nab,1);

  t_norm=my_dnrm2(ND,t_nbc,1);
  if (t_norm < D_EPS) {
    *omega=0.0;
    return(0);
  }
  my_dscal(ND,1.0/t_norm,t_nbc,1);

  t_norm=my_dnrm2(ND,t_nca,1);
  if (t_norm < D_EPS) {
    *omega=0.0;
    return(0);
  }
  my_dscal(ND,1.0/t_norm,t_nca,1);

  /* calculate dihedral angles between facets */
  /* TODO source of this formula ? */

  PetscReal t_a_abbc,t_a_bcca,t_a_caab;
  t_a_abbc=t_nab[0]*t_nbc[0]+t_nab[1]*t_nbc[1]+t_nab[2]*t_nbc[2];
  t_a_bcca=t_nbc[0]*t_nca[0]+t_nbc[1]*t_nca[1]+t_nbc[2]*t_nca[2];
  t_a_caab=t_nca[0]*t_nab[0]+t_nca[1]*t_nab[1]+t_nca[2]*t_nab[2];
  if (t_a_abbc>1) t_a_abbc=PETSC_PI; else if (t_a_abbc<-1) t_a_abbc=0; else t_a_abbc=PETSC_PI-acos(t_nab[0]*t_nbc[0]+t_nab[1]*t_nbc[1]+t_nab[2]*t_nbc[2]);
  if (t_a_bcca>1) t_a_bcca=PETSC_PI; else if (t_a_bcca<-1) t_a_bcca=0; else t_a_bcca=PETSC_PI-acos(t_nbc[0]*t_nca[0]+t_nbc[1]*t_nca[1]+t_nbc[2]*t_nca[2]);
  if (t_a_caab>1) t_a_caab=PETSC_PI; else if (t_a_caab<-1) t_a_caab=0; else t_a_caab=PETSC_PI-acos(t_nca[0]*t_nab[0]+t_nca[1]*t_nab[1]+t_nca[2]*t_nab[2]);

  *omega=t_a_abbc+t_a_bcca+t_a_caab-PETSC_PI;

  return(0);
}