Skip to content

Fix a critical bug in scuff-tmatrix (undesired VSWF scaling) #197

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
texnokrates wants to merge 8 commits into
base: master
Choose a base branch
from
Open

Fix a critical bug in scuff-tmatrix (undesired VSWF scaling) #197

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5f64c76
add fclose to RWGGeometry constructor
Jan 20, 2019
69e9411
const hygiene in libSpherical coordinate transforms.
texnokrates Apr 5, 2019
97751ba
Prune some effectively no-op code.
texnokrates Apr 5, 2019
6403bdf
libSpherical.cc: Remove auxiliary variables to improve readability.
texnokrates May 21, 2019
78689f5
Fix scuff-tmatrix: Add missing parameter, preventing unwanted VSW sca…
texnokrates May 23, 2019
32f080f
Merge branch 'addfclose' of https://github.com/WenjieYao/scuff-em
texnokrates May 24, 2019
cb1109b
Determine version from git
texnokrates May 24, 2019
8e7562c
scuff-tmatrix: write version to output
texnokrates May 24, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions applications/scuff-tmatrix/Makefile.am
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,10 @@ scuff_tmatrix_SOURCES = \

scuff_tmatrix_LDADD = $(top_builddir)/libs/libscuff/libscuff.la

VERSIONSTRING=\"$(shell git describe --dirty)\"

AM_CPPFLAGS = -DSCUFF \
-DVERSION=$(VERSIONSTRING) \
-I$(top_srcdir)/libs/libscuff \
-I$(top_srcdir)/libs/libIncField \
-I$(top_srcdir)/libs/libMatProp \
Expand Down
8 changes: 7 additions & 1 deletion applications/scuff-tmatrix/scuff-tmatrix.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,12 @@
* http://homerreid.github.io/scuff-em-documentation/applications/scuff-tmatrix/scuff-tmatrix
*/

#ifdef VERSION
#define VERSIONSTRING VERSION " "
#else
#define VERSIONSTRING ""
#endif

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

Expand Down Expand Up @@ -133,7 +139,7 @@ int main(int argc, char *argv[])
if (!FileBase)
FileBase = vstrdup(GetFileBase(GeoFileName));
FILE *f=vfopen("%s.TMatrix","a",FileBase);
fprintf(f,"# scuff-tmatrix run on %s (%s)\n",GetHostName(),GetTimeString());
fprintf(f,"# scuff-tmatrix " VERSIONSTRING "run on %s (%s)\n",GetHostName(),GetTimeString());
fprintf(f,"# columns:\n");
fprintf(f,"# 1 omega\n");
fprintf(f,"# 2,3,4,5 (alpha, {L,M,P}_alpha) (T-matrix row index)\n");
Expand Down
2 changes: 1 addition & 1 deletion configure.ac
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
AC_INIT(scuff-em, 0.96, homer@homerreid.com)
AC_INIT(scuff-em, [m4_esyscmd_s([git describe --always --dirty])], homer@homerreid.com)
AC_CONFIG_SRCDIR([README])

# Shared-library version. This has a technical
Expand Down
7 changes: 1 addition & 6 deletions libs/libIncField/SphericalWave.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,15 +23,12 @@
* homer reid -- 11/2009 -- 2/2012
*/

#include <string.h>
#include <math.h>

#include <libhrutil.h>
#include <libIncField.h>
#include <libSpherical.h>

#define II cdouble(0.0,1.0)

/**********************************************************************/
/* constructor and field-setting routines *****************************/
/**********************************************************************/
Expand All @@ -57,10 +54,8 @@ void SphericalWave::GetFields(const double X[3], cdouble EHC[6])
cdouble Z=ZVAC*sqrt(Mu/Eps);

// convert the evaluation point to spherical coordinates
double XX[3];
memcpy(XX, X, 3*sizeof(double));
double r, Theta, Phi;
CoordinateC2S(XX, &r, &Theta, &Phi);
CoordinateC2S(X, &r, &Theta, &Phi);

// get the M and N vector spherical harmonics
cdouble MVec[3], NVec[3];
Expand Down
2 changes: 1 addition & 1 deletion libs/libSpherical/VectorSphericalWaves.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -246,7 +246,7 @@ HMatrix *GetWaveMatrix(double r, double Theta, double Phi,
cdouble *RFArray = dYdThetaArray + NAlpha;
cdouble *dWorkspace = RFArray + 3*(LMax+1);

GetVSWRadialFunctions(LMax, k, r, WaveType, RFArray, (double *)dWorkspace, RConjugate);
GetVSWRadialFunctions(LMax, k, r, WaveType, RFArray, (double *)dWorkspace, /* TimesrFactor */ false, RConjugate);

/***************************************************************/
/* fetch angular functions */
Expand Down
33 changes: 15 additions & 18 deletions libs/libSpherical/libSpherical.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,10 +52,10 @@ void CoordinateC2S(double X, double Y, double Z, double *r, double *Theta, doubl
*Phi=atan2(Y,X);
}

void CoordinateC2S(double X[3], double *r, double *Theta, double *Phi)
void CoordinateC2S(const double X[3], double *r, double *Theta, double *Phi)
{ CoordinateC2S(X[0], X[1], X[2], r, Theta, Phi); }

void CoordinateC2S(double X[3], double R[3])
void CoordinateC2S(const double X[3], double R[3])
{ CoordinateC2S(X[0], X[1], X[2], R+0, R+1, R+2); }

void CoordinateC2S(double X[3])
Expand All @@ -82,7 +82,7 @@ void CoordinateS2C(double r, double Theta, double Phi, double *X, double *Y, dou
void CoordinateS2C(double r, double Theta, double Phi, double X[3])
{ CoordinateS2C(r, Theta, Phi, X+0, X+1, X+2); }

void CoordinateS2C(double R[3], double X[3])
void CoordinateS2C(const double R[3], double X[3])
{ CoordinateS2C(R[0], R[1], R[2], X+0, X+1, X+2); }

void CoordinateS2C(double R[3])
Expand All @@ -97,7 +97,7 @@ void CoordinateS2C(double R[3])
/* given the cartesian components of a vector, return its */
/* spherical components. */
/***************************************************************/
void VectorC2S(double Theta, double Phi, cdouble VC[3], cdouble VS[3])
void VectorC2S(double Theta, double Phi, const cdouble VC[3], cdouble VS[3])
{
double CT=cos(Theta), ST=sin(Theta);
double CP=cos(Phi), SP=sin(Phi);
Expand Down Expand Up @@ -126,7 +126,7 @@ void VectorC2S(double Theta, double Phi, double V[3])
V[0]=VS[0]; V[1]=VS[1]; V[2]=VS[2];
}

void VectorC2S(double Theta, double Phi, double VC[3], double VS[3])
void VectorC2S(double Theta, double Phi, const double VC[3], double VS[3])
{
double CT=cos(Theta), ST=sin(Theta);
double CP=cos(Phi), SP=sin(Phi);
Expand Down Expand Up @@ -160,7 +160,7 @@ void VectorS2C(double Theta, double Phi, double V[3])
/* given the spherical components of a vector, return its */
/* cartesian components. */
/***************************************************************/
void VectorS2C(double Theta, double Phi, cdouble VS[3], cdouble VC[3])
void VectorS2C(double Theta, double Phi, const cdouble VS[3], cdouble VC[3])
{
double CT=cos(Theta), ST=sin(Theta);
double CP=cos(Phi), SP=sin(Phi);
Expand All @@ -178,7 +178,7 @@ void VectorS2C(double Theta, double Phi, cdouble VS[3], cdouble VC[3])

}

void VectorS2C(double Theta, double Phi, double VS[3], double VC[3])
void VectorS2C(double Theta, double Phi, const double VS[3], double VC[3])
{
double CT=cos(Theta), ST=sin(Theta);
double CP=cos(Phi), SP=sin(Phi);
Expand Down Expand Up @@ -760,7 +760,7 @@ void GetMNlmArray(int lMax, cdouble k,
{
int NAlpha=(lMax+1)*(lMax+1);
int Alpha, l, m;
double SinTheta, dl, dm;
double SinTheta;
cdouble nik, PreFac;
cdouble ROverR;
cdouble *R = new cdouble[lMax+2];
Expand Down Expand Up @@ -798,13 +798,10 @@ void GetMNlmArray(int lMax, cdouble k,
for (Alpha=l=1; l<=lMax; l++)
for (m=-l; m<=l; m++, Alpha++)
{
dl=((double)l);
dm=((double)m);

PreFac=1.0/sqrt( dl*(dl+1.0) );
PreFac=1.0/sqrt( l*(l+1.0) );

M[3*Alpha + 0]= 0.0;
M[3*Alpha + 1]= -dm*PreFac*R[l]*Ylm[Alpha]/SinTheta;
M[3*Alpha + 1]= -(double)m*PreFac*R[l]*Ylm[Alpha]/SinTheta;
M[3*Alpha + 2]= -II*PreFac*R[l]*dYlmdTheta[Alpha];

#if 0
Expand All @@ -815,24 +812,24 @@ void GetMNlmArray(int lMax, cdouble k,
ROverR = (l==1) ? k/3.0 : 0.0;
else
ROverR=R[l]/r;
N[3*Alpha + 0]= -sqrt(dl*(dl+1.0))*ROverR*Ylm[Alpha]/k;
N[3*Alpha + 0]= -sqrt(l*(l+1.0))*ROverR*Ylm[Alpha]/k;
N[3*Alpha + 1]= II*PreFac*(ROverR + dRdr[l])*dYlmdTheta[Alpha];
N[3*Alpha + 2]= -dm*PreFac*(ROverR + dRdr[l])*Ylm[Alpha]/SinTheta;
N[3*Alpha + 2]= -(double)m*PreFac*(ROverR + dRdr[l])*Ylm[Alpha]/SinTheta;
#endif
PreFac/=k;

if (r==0.0)
ROverR = (l==1) ? k/3.0 : 0.0;
else
ROverR=R[l]/r;
N[3*Alpha + 0]= -sqrt(dl*(dl+1.0))*ROverR*Ylm[Alpha]/k;
N[3*Alpha + 0]= -sqrt(l*(l+1.0))*ROverR*Ylm[Alpha]/k;
N[3*Alpha + 1]= -PreFac*(ROverR + dRdr[l])*dYlmdTheta[Alpha];
N[3*Alpha + 2]= -II*dm*PreFac*(ROverR + dRdr[l])*Ylm[Alpha]/SinTheta;
N[3*Alpha + 2]= -II*(double)m*PreFac*(ROverR + dRdr[l])*Ylm[Alpha]/SinTheta;

if (LL)
{ LL[3*Alpha + 0] = PreFac*dRdr[l] * Ylm[Alpha];
LL[3*Alpha + 1] = PreFac*ROverR*dYlmdTheta[Alpha];
LL[3*Alpha + 2] = PreFac*II*dm*ROverR * Ylm[Alpha] / SinTheta;
LL[3*Alpha + 2] = PreFac*II*(double)m*ROverR * Ylm[Alpha] / SinTheta;
};

if (DivLL)
Expand Down
16 changes: 8 additions & 8 deletions libs/libSpherical/libSpherical.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,19 +67,19 @@
/* conversion routines */
/***************************************************************/
void CoordinateC2S(double X, double Y, double Z, double *r, double *Theta, double *Phi);
void CoordinateC2S(double X[3], double *r, double *Theta, double *Phi);
void CoordinateC2S(double X[3], double R[3]);
void CoordinateC2S(double X[3]);
void CoordinateC2S(const double X[3], double *r, double *Theta, double *Phi);
void CoordinateC2S(const double X[3], double R[3]);
void CoordinateC2S(const double X[3]);
void CoordinateS2C(double r, double Theta, double Phi, double X[3]);
void CoordinateS2C(double r, double Theta, double Phi, double *X, double *Y, double *Z);
void CoordinateS2C(double R[3], double X[3]);
void CoordinateS2C(const double R[3], double X[3]);
void CoordinateS2C(double R[3]);
void VectorC2S(double Theta, double Phi, cdouble VC[3], cdouble VS[3]);
void VectorC2S(double Theta, double Phi, double VC[3], double VS[3]);
void VectorC2S(double Theta, double Phi, const cdouble VC[3], cdouble VS[3]);
void VectorC2S(double Theta, double Phi, const double VC[3], double VS[3]);
void VectorC2S(double Theta, double Phi, double V[3]);
void VectorC2S(double Theta, double Phi, cdouble V[3]);
void VectorS2C(double Theta, double Phi, cdouble VS[3], cdouble VC[3]);
void VectorS2C(double Theta, double Phi, double VS[3], double VC[3]);
void VectorS2C(double Theta, double Phi, const cdouble VS[3], cdouble VC[3]);
void VectorS2C(double Theta, double Phi, const double VS[3], double VC[3]);
void VectorS2C(double Theta, double Phi, double V[3]);
void VectorS2C(double Theta, double Phi, cdouble V[3]);

Expand Down
1 change: 1 addition & 0 deletions libs/libscuff/RWGGeometry.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -804,6 +804,7 @@ RWGGeometry::RWGGeometry(const char *pGeoFileName, int pLogLevel)
FIBBICaches[ns] = FIBBICaches[ Mate[ns] ];
else
FIBBICaches[ns] = CreateFIBBICache(Surfaces[ns]->MeshFileName);
fclose(f);
}

/***************************************************************/
Expand Down