util_sparx.cpp File Reference

#include <cstring>
#include <ctime>
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <boost/format.hpp>
#include "emdata.h"
#include "util.h"
#include "fundamentals.h"
#include "lapackblas.h"
#include "lbfgsb.h"
#include "steepest.h"
#include "emassert.h"
#include "randnum.h"
#include <gsl/gsl_sf_bessel.h>
#include <cmath>

Include dependency graph for util_sparx.cpp:

Go to the source code of this file.

Classes
struct	ori_t
struct	cmpang
struct	tmpstruct
struct	stcom_
struct	ccf_point
struct	ccf_value
struct	point3d_t
Defines
#define	fdata(i, j)   fdata[ i-1 + (j-1)*nxdata ]
#define	fdata(i, j)   fdata[ i-1 + (j-1)*nxdata ]
#define	circ(i)   circ[i-1]
#define	numr(i, j)   numr[(j-1)*3 + i-1]
#define	xim(i, j)   xim[(j-1)*nsam + i-1]
#define	tab1(i)   tab1[i-1]
#define	xcmplx(i, j)   xcmplx [(j-1)*2 + i-1]
#define	br(i)   br[i-1]
#define	bi(i)   bi[i-1]
#define	b(i)   b[i-1]
#define	circ1(i)   circ1[i-1]
#define	circ2(i)   circ2[i-1]
#define	t(i)   t[i-1]
#define	q(i)   q[i-1]
#define	b(i)   b[i-1]
#define	t7(i)   t7[i-1]
#define	dout(i, j)   dout[i+maxrin*j]
#define	circ1b(i)   circ1b[i-1]
#define	circ2b(i)   circ2b[i-1]
#define	dout(i, j)   dout[i+maxrin*j]
#define	circ1b(i)   circ1b[i-1]
#define	circ2b(i)   circ2b[i-1]
#define	QUADPI   3.141592653589793238462643383279502884197
#define	PI2   2*QUADPI
#define	QUADPI   3.141592653589793238462643383279502884197
#define	PI2   QUADPI*2
#define	deg_rad   QUADPI/180.0
#define	rad_deg   180.0/QUADPI
#define	old_ptr(i, j, k)   old_ptr[i+(j+(k*ny))*nx]
#define	new_ptr(iptr, jptr, kptr)   new_ptr[iptr+(jptr+(kptr*new_ny))*new_nx]
#define	inp(i, j, k)   inp[(i+new_st_x)+((j+new_st_y)+((k+new_st_z)*ny))*nx]
#define	outp(i, j, k)   outp[i+(j+(k*new_ny))*new_nx]
#define	inp(i, j, k)   inp[i+(j+(k*ny))*nx]
#define	outp(i, j, k)   outp[(i+new_st_x)+((j+new_st_y)+((k+new_st_z)*new_ny))*new_nx]
#define	QUADPI   3.141592653589793238462643383279502884197
#define	DGR_TO_RAD   QUADPI/180
#define	DM(I)   DM [I-1]
#define	SS(I)   SS [I-1]
#define	DM(I)   DM[I-1]
#define	B(i, j)   Bptr[i-1+((j-1)*NSAM)]
#define	CUBE(i, j, k)   CUBEptr[(i-1)+((j-1)+((k-1)*NY3D))*NX3D]
#define	W(i, j)   Wptr [i-1+((j-1)*Wnx)]
#define	PROJ(i, j)   PROJptr [i-1+((j-1)*NNNN)]
#define	SS(I, J)   SS [I-1 + (J-1)*6]
#define	W(i, j)   Wptr [i-1+((j-1)*Wnx)]
#define	PROJ(i, j)   PROJptr [i-1+((j-1)*NNNN)]
#define	SS(I, J)   SS [I-1 + (J-1)*6]
#define	RI(i, j)   RI [(i-1) + ((j-1)*3)]
#define	CC(i)   CC [i-1]
#define	CP(i)   CP [i-1]
#define	VP(i)   VP [i-1]
#define	VV(i)   VV [i-1]
#define	AMAX1(i, j)   i>j?i:j
#define	AMIN1(i, j)   i<j?i:j
#define	mymax(x, y)   (((x)>(y))?(x):(y))
#define	mymin(x, y)   (((x)<(y))?(x):(y))
#define	sign(x, y)   (((((y)>0)?(1):(-1))*(y!=0))*(x))
#define	quadpi   3.141592653589793238462643383279502884197
#define	dgr_to_rad   quadpi/180
#define	deg_to_rad   quadpi/180
#define	rad_to_deg   180/quadpi
#define	rad_to_dgr   180/quadpi
#define	TRUE   1
#define	FALSE   0
#define	theta(i)   theta [i-1]
#define	phi(i)   phi [i-1]
#define	weight(i)   weight [i-1]
#define	lband(i)   lband [i-1]
#define	ts(i)   ts [i-1]
#define	thetast(i)   thetast [i-1]
#define	key(i)   key [i-1]
#define	TRUE_   (1)
#define	FALSE_   (0)
#define	abs(x)   ((x) >= 0 ? (x) : -(x))
#define	img_ptr(i, j, k)   img_ptr[2*(i-1)+((j-1)+((k-1)*ny))*nxo]
#define	img_ptr(i, j, k)   img_ptr[i+(j+(k*ny))*nx]
#define	img2_ptr(i, j, k)   img2_ptr[i+(j+(k*ny))*nx]
#define	cent(i)   out[i+N]
#define	assign(i)   out[i]
#define	data(i, j)   group[i*ny+j]
Functions
int	i_dnnt (double *x)
int	addnod_ (int *nst, int *k, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *lnew, int *ier)
double	angle_ (double *v1, double *v2, double *v3)
double	areas_ (double *v1, double *v2, double *v3)
double	areav_new__ (int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *ier)
int	bdyadd_ (int *kk, int *i1, int *i2, int *list, int *lptr, int *lend, int *lnew)
int	bnodes_ (int *n, int *list, int *lptr, int *lend, int *nodes, int *nb, int *na, int *nt)
int	circle_ (int *k, double *xc, double *yc, int *ier)
int	circum_ (double *v1, double *v2, double *v3, double *c__, int *ier)
int	covsph_ (int *kk, int *n0, int *list, int *lptr, int *lend, int *lnew)
int	crlist_ (int *n, int *ncol, double *x, double *y, double *z__, int *list, int *lend, int *lptr, int *lnew, int *ltri, int *listc, int *nb, double *xc, double *yc, double *zc, double *rc, int *ier)
int	delarc_ (int *n, int *io1, int *io2, int *list, int *lptr, int *lend, int *lnew, int *ier)
int	delnb_ (int *n0, int *nb, int *n, int *list, int *lptr, int *lend, int *lnew, int *lph)
int	delnod_ (int *k, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *lnew, int *lwk, int *iwk, int *ier)
int	drwarc_ (int *, double *p, double *q, double *tol, int *nseg)
int	edge_ (int *in1, int *in2, double *x, double *y, double *z__, int *lwk, int *iwk, int *list, int *lptr, int *lend, int *ier)
int	getnp_ (double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *l, int *npts, double *df, int *ier)
int	insert_ (int *k, int *lp, int *list, int *lptr, int *lnew)
long int	inside_ (double *p, int *lv, double *xv, double *yv, double *zv, int *nv, int *listv, int *ier)
int	intadd_ (int *kk, int *i1, int *i2, int *i3, int *list, int *lptr, int *lend, int *lnew)
int	intrsc_ (double *p1, double *p2, double *cn, double *p, int *ier)
int	jrand_ (int *n, int *ix, int *iy, int *iz)
long int	left_ (double *x1, double *y1, double *z1, double *x2, double *y2, double *z2, double *x0, double *y0, double *z0)
int	lstptr_ (int *lpl, int *nb, int *list, int *lptr)
int	nbcnt_ (int *lpl, int *lptr)
int	nearnd_ (double *p, int *ist, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, double *al)
int	optim_ (double *x, double *y, double *z__, int *na, int *list, int *lptr, int *lend, int *nit, int *iwk, int *ier)
int	projct_ (double *px, double *py, double *pz, double *ox, double *oy, double *oz, double *ex, double *ey, double *ez, double *vx, double *vy, double *vz, long int *init, double *x, double *y, double *z__, int *ier)
int	scoord_ (double *px, double *py, double *pz, double *plat, double *plon, double *pnrm)
double	store_ (double *x)
int	swap_ (int *in1, int *in2, int *io1, int *io2, int *list, int *lptr, int *lend, int *lp21)
long int	swptst_ (int *n1, int *n2, int *n3, int *n4, double *x, double *y, double *z__)
int	trans_ (int *n, double *rlat, double *rlon, double *x, double *y, double *z__)
int	trfind_ (int *nst, double *p, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, double *b1, double *b2, double *b3, int *i1, int *i2, int *i3)
int	trlist_ (int *n, int *list, int *lptr, int *lend, int *nrow, int *nt, int *ltri, int *ier)
int	trlprt_ (int *n, double *x, double *y, double *z__, int *iflag, int *nrow, int *nt, int *ltri, int *lout)
int	trmesh_ (int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, int *lnew, int *near__, int *next, double *dist, int *ier)
int	trplot_ (int *lun, double *pltsiz, double *elat, double *elon, double *a, int *n, double *x, double *y, double *z__, int *list, int *lptr, int *lend, char *, long int *numbr, int *ier, short)
int	trprnt_ (int *n, double *x, double *y, double *z__, int *iflag, int *list, int *lptr, int *lend, int *lout)
int	vrplot_ (int *lun, double *pltsiz, double *elat, double *elon, double *a, int *n, double *x, double *y, double *z__, int *nt, int *listc, int *lptr, int *lend, double *xc, double *yc, double *zc, char *, long int *numbr, int *ier, short)
int	random_ (int *ix, int *iy, int *iz, double *rannum)
int	find_group (int ix, int iy, int iz, int grpid, EMData *mg, EMData *visited)
Variables
stcom_	stcom_1

---

Define Documentation

#define abs

(

x

)

((x) >= 0 ? (x) : -(x))

Definition at line 7502 of file util_sparx.cpp.

#define AMAX1	(	i,
		j		)	i>j?i:j

Definition at line 5617 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define AMIN1	(	i,
		j		)	i<j?i:j

Definition at line 5618 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define assign

(

i

)

out[i]

Definition at line 18926 of file util_sparx.cpp.

Referenced by EMAN::Util::cluster_pairwise().

#define B	(	i,
		j		)	Bptr[i-1+((j-1)*NSAM)]

Definition at line 5463 of file util_sparx.cpp.

Referenced by EMAN::Util::BPCQ(), EMAN::Util::histc(), EMAN::Util::im_diff(), and submatrix().

#define b

(

i

)

b[i-1]

Definition at line 3189 of file util_sparx.cpp.

#define b

(

i

)

b[i-1]

Definition at line 3189 of file util_sparx.cpp.

Referenced by EMAN::CtfCWautoAverager::add_image(), bmv_(), EMAN::Util::cml_line_insino(), EMAN::Util::cml_line_insino_all(), EMAN::OptVarianceCmp::cmp(), Derivatives(), Derivatives_G(), dpmeps_(), formk_(), GCVmin_Tik(), EMAN::TetrahedralSym::get_asym_unit_points(), EMAN::PlatonicSym::get_asym_unit_points(), EMAN::HSym::get_asym_unit_points(), EMAN::EMUtil::get_euler_names(), inside_(), EMAN::Matrix4::inverse(), main(), EMAN::Matrix4::operator*(), EMAN::Quaternion::operator*=(), EMAN::Quaternion::operator/=(), EMAN::Util::prb1d(), prb1d(), EMAN::TestImageEllipse::process_inplace(), EMAN::TestImageGradient::process_inplace(), EMAN::NormalizeToLeastSquareProcessor::process_inplace(), EMAN::GradientRemoverProcessor::process_inplace(), EMAN::Util::splint(), subsm_(), and varmx().

#define bi

(

i

)

bi[i-1]

Definition at line 2642 of file util_sparx.cpp.

Referenced by EMAN::Util::fftc_d(), fftc_d(), EMAN::Util::fftc_q(), fftc_q(), and EMAN::TestImageEllipse::process_inplace().

#define br

(

i

)

br[i-1]

Definition at line 2641 of file util_sparx.cpp.

Referenced by EMAN::Util::fftc_d(), fftc_d(), EMAN::Util::fftc_q(), fftc_q(), EMAN::EMData::render_amp24(), and EMAN::EMData::render_ap24().

#define CC

(

i

)

CC [i-1]

Definition at line 5613 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define cent

(

i

)

out[i+N]

Definition at line 18925 of file util_sparx.cpp.

Referenced by EMAN::Util::cluster_pairwise().

#define circ

(

i

)

circ[i-1]

Definition at line 2159 of file util_sparx.cpp.

Referenced by EMAN::Util::alrl_ms(), alrq(), alrq_ms(), applyws(), EMAN::Util::Frngs(), frngs(), EMAN::Util::Frngs_inv(), EMAN::Util::Polar2D(), EMAN::Util::Polar2Dm(), and EMAN::Util::Polar2Dmi().

#define circ1

(

i

)

circ1[i-1]

Definition at line 3185 of file util_sparx.cpp.

Referenced by EMAN::Util::Crosrng_e(), crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), crosrng_ms(), EMAN::Util::Crosrng_ns(), and EMAN::Util::Crosrng_sm_psi().

#define circ1b

(

i

)

circ1b[i-1]

Definition at line 3912 of file util_sparx.cpp.

#define circ1b

(

i

)

circ1b[i-1]

Definition at line 3912 of file util_sparx.cpp.

Referenced by EMAN::Util::Crosrng_msg(), EMAN::Util::Crosrng_msg_m(), EMAN::Util::Crosrng_msg_s(), and EMAN::Util::Crosrng_msg_vec().

#define circ2

(

i

)

circ2[i-1]

Definition at line 3186 of file util_sparx.cpp.

Referenced by EMAN::Util::Crosrng_e(), crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), crosrng_ms(), EMAN::Util::Crosrng_ns(), and EMAN::Util::Crosrng_sm_psi().

#define circ2b

(

i

)

circ2b[i-1]

Definition at line 3913 of file util_sparx.cpp.

#define circ2b

(

i

)

circ2b[i-1]

Definition at line 3913 of file util_sparx.cpp.

Referenced by EMAN::Util::Crosrng_msg(), EMAN::Util::Crosrng_msg_m(), EMAN::Util::Crosrng_msg_s(), and EMAN::Util::Crosrng_msg_vec().

#define CP

(

i

)

CP [i-1]

Definition at line 5614 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define CUBE	(	i,
		j,
		k		)	CUBEptr[(i-1)+((j-1)+((k-1)*NY3D))*NX3D]

Definition at line 5464 of file util_sparx.cpp.

Referenced by EMAN::Util::BPCQ().

#define data	(	i,
		j		)	group[i*ny+j]

Definition at line 19232 of file util_sparx.cpp.

Referenced by EMAN::EMData::absi(), EMAN::EMData::add(), EMAN::file_store::add_image(), EMAN::EMData::addsquare(), EMAN::RotatePrecenterAligner::align(), EMAN::RotationalAligner::align_180_ambiguous(), EMAN::EMData::amplitude(), EMAN::EMData::apply_radial_func(), EMAN::EMData::calc_az_dist(), EMAN::EMData::calc_center_of_mass(), EMAN::EMData::calc_highest_locations(), EMAN::EMData::calc_hist(), EMAN::MaskEdgeMeanProcessor::calc_locals(), EMAN::EMData::calc_max_location(), EMAN::EMData::calc_min_location(), EMAN::EMData::calc_radial_dist(), circumference(), EMAN::BoxingTools::classify(), EMAN::EMData::common_lines(), EMAN::EMData::common_lines_real(), EMAN::Util::cyclicshift(), EMAN::PointArray::distmx(), EMAN::EMData::div(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::EMData(), EMAN::EMData::get_attr(), EMAN::EMData::get_circle_mean(), get_data_as_vector(), EMAN::EMData::get_edge_mean(), EMAN::EMData::get_fft_amplitude(), EMAN::EMData::get_fft_phase(), EMAN::file_store::get_image(), EMAN::newfile_store::get_image(), EMAN::Util::histc(), EMAN::EMData::imag(), EMAN::EMData::insert_scaled_sum(), EMAN::SingleSpiderIO::is_valid(), EMAN::SpiderIO::is_valid(), EMAN::PifIO::is_valid(), EMAN::MrcIO::is_valid(), EMAN::ImagicIO::is_valid(), EMAN::IcosIO::is_valid(), EMAN::Gatan2IO::is_valid(), EMAN::EmIO::is_valid(), EMAN::EmimIO::is_valid(), EMAN::DM3IO::is_valid(), EMAN::EMData::little_big_dot(), EMAN::EMData::log(), EMAN::EMData::log10(), main(), EMAN::TestUtil::make_image_file_by_mode(), mpi_init(), mpi_start(), EMAN::EMData::mult(), EMAN::EMData::mult_complex_efficient(), EMAN::EMData::norm_pad(), EMAN::Util::Normalize_ring(), EMAN::EMData::operator=(), EMAN::PointArray::pdb2mrc_by_nfft(), EMAN::EMData::phase(), EMAN::XYZProcessor::process_inplace(), EMAN::CutoffBlockProcessor::process_inplace(), EMAN::DiffBlockProcessor::process_inplace(), EMAN::BoxStatProcessor::process_inplace(), EMAN::AreaProcessor::process_inplace(), EMAN::ComplexPixelProcessor::process_inplace(), EMAN::CoordinateProcessor::process_inplace(), EMAN::RealPixelProcessor::process_inplace(), EMAN::ImageProcessor::process_inplace(), EMAN::BoxMedianProcessor::process_pixel(), EMAN::GaussFFTProjector::project3d(), EMAN::PointArray::projection_by_nfft(), EMAN::Gatan::TagData::read_array_data(), EMAN::EMData::real(), EMAN::EMData::render_amp24(), EMAN::EMData::render_ap24(), EMAN::EMData::ri2ap(), EMAN::EMData::ri2inten(), EMAN::EMData::rot_scale_conv_new(), EMAN::EMData::rot_scale_conv_new_background(), EMAN::EMData::rotate_x(), EMAN::BoxSVDClassifier::setDims(), EMAN::EMData::setup4slice(), EMAN::EMData::sqrt(), EMAN::EMData::sub(), EMAN::EMData::subsquare(), EMAN::EMData::to_value(), EMAN::EMData::update_stat(), EMAN::Util::vareas(), EMAN::TestUtil::verify_image_file_by_mode(), EMAN::EMUtil::vertical_acf(), and EMAN::U3DWriter::write_clod_mesh_generator_node().

#define deg_rad   QUADPI/180.0

Definition at line 4322 of file util_sparx.cpp.

Referenced by EMAN::Util::cml_init_rot(), EMAN::Util::cml_line_in3d(), and EMAN::Util::cml_update_rot().

#define deg_to_rad   quadpi/180

Definition at line 6793 of file util_sparx.cpp.

#define dgr_to_rad   quadpi/180

Definition at line 6792 of file util_sparx.cpp.

Referenced by EMAN::Util::ang_to_xyz(), apmq(), aprq2d(), EMAN::Util::even_angles(), and EMAN::ChaoProjector::setdm().

#define DGR_TO_RAD   QUADPI/180

Definition at line 5415 of file util_sparx.cpp.

#define DM

(

I

)

DM[I-1]

Definition at line 5462 of file util_sparx.cpp.

#define DM

(

I

)

DM [I-1]

Definition at line 5462 of file util_sparx.cpp.

Referenced by EMAN::Util::BPCQ(), and EMAN::Util::CANG().

#define dout	(	i,
		j		)	dout[i+maxrin*j]

Definition at line 3911 of file util_sparx.cpp.

#define dout	(	i,
		j		)	dout[i+maxrin*j]

Definition at line 3911 of file util_sparx.cpp.

Referenced by EMAN::Util::Crosrng_msg(), EMAN::Util::Crosrng_msg_m(), and EMAN::Util::Crosrng_msg_s().

#define FALSE   0

Definition at line 6797 of file util_sparx.cpp.

#define FALSE_   (0)

Definition at line 7501 of file util_sparx.cpp.

#define fdata	(	i,
		j		)	fdata[ i-1 + (j-1)*nxdata ]

Definition at line 709 of file util_sparx.cpp.

#define fdata	(	i,
		j		)	fdata[ i-1 + (j-1)*nxdata ]

Definition at line 709 of file util_sparx.cpp.

Referenced by EMAN::Util::quadri(), quadri(), and EMAN::Util::quadri_background().

#define img2_ptr	(	i,
		j,
		k		)	img2_ptr[i+(j+(k*ny))*nx]

Definition at line 18550 of file util_sparx.cpp.

Referenced by EMAN::Util::addn_img(), EMAN::Util::divn_filter(), EMAN::Util::divn_img(), EMAN::Util::madn_scalar(), EMAN::Util::move_points(), EMAN::Util::muln_img(), EMAN::Util::mult_scalar(), and EMAN::Util::subn_img().

#define img_ptr	(	i,
		j,
		k		)	img_ptr[i+(j+(k*ny))*nx]

Definition at line 18549 of file util_sparx.cpp.

#define img_ptr	(	i,
		j,
		k		)	img_ptr[2*(i-1)+((j-1)+((k-1)*ny))*nxo]

Definition at line 18549 of file util_sparx.cpp.

Referenced by EMAN::Util::add_img(), EMAN::Util::add_img2(), EMAN::Util::addn_img(), EMAN::Util::compress_image_mask(), EMAN::Util::div_filter(), EMAN::Util::div_img(), EMAN::Util::divn_filter(), EMAN::Util::divn_img(), EMAN::Util::hist_comp_freq(), EMAN::Util::mad_scalar(), EMAN::Util::madn_scalar(), EMAN::Util::move_points(), EMAN::Util::mul_img(), EMAN::Util::mul_scalar(), EMAN::Util::muln_img(), EMAN::Util::mult_scalar(), EMAN::Util::pack_complex_to_real(), ReadStackandDist(), ReadStackandDist_Cart(), EMAN::Util::reconstitute_image_mask(), EMAN::Util::set_line(), EMAN::Util::sub_img(), and EMAN::Util::subn_img().

#define inp	(	i,
		j,
		k		)	inp[i+(j+(k*ny))*nx]

Definition at line 5049 of file util_sparx.cpp.

#define inp	(	i,
		j,
		k		)	inp[(i+new_st_x)+((j+new_st_y)+((k+new_st_z)*ny))*nx]

Definition at line 5049 of file util_sparx.cpp.

Referenced by EMAN::Util::pad(), and EMAN::Util::window().

#define key

(

i

)

key [i-1]

Definition at line 6806 of file util_sparx.cpp.

Referenced by EMAN::Dict::get(), EMAN::Util::hsortd(), mpi_comm_split(), EMAN::Log::vlog(), EMAN::Util::voronoi(), and EMAN::Util::vrdg().

#define lband

(

i

)

lband [i-1]

Definition at line 6803 of file util_sparx.cpp.

#define mymax	(	x,
		y		)	(((x)>(y))?(x):(y))

Definition at line 6786 of file util_sparx.cpp.

#define mymin	(	x,
		y		)	(((x)<(y))?(x):(y))

Definition at line 6787 of file util_sparx.cpp.

#define new_ptr	(	iptr,
		jptr,
		kptr		)	new_ptr[iptr+(jptr+(kptr*new_ny))*new_nx]

Definition at line 4945 of file util_sparx.cpp.

Referenced by EMAN::Util::compress_image_mask(), EMAN::Util::decimate(), and EMAN::Util::reconstitute_image_mask().

#define numr	(	i,
		j		)	numr[(j-1)*3 + i-1]

Definition at line 2160 of file util_sparx.cpp.

Referenced by ali3d_d(), alprbs(), EMAN::Util::alrl_ms(), alrq(), alrq_ms(), apmd(), apmq(), applyws(), apring1(), aprq2d(), EMAN::Util::Crosrng_e(), crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), crosrng_ms(), EMAN::Util::Crosrng_msg(), EMAN::Util::Crosrng_msg_m(), EMAN::Util::Crosrng_msg_s(), EMAN::Util::Crosrng_msg_vec(), EMAN::Util::Crosrng_ns(), EMAN::Util::Crosrng_sm_psi(), EMAN::Util::ener(), EMAN::Util::ener_tot(), EMAN::Util::Frngs(), frngs(), EMAN::Util::Frngs_inv(), numrinit(), Numrinit(), EMAN::Util::Polar2D(), EMAN::Util::Polar2Dm(), EMAN::Util::Polar2Dmi(), ringwe(), EMAN::Util::sub_fav(), and EMAN::Util::update_fav().

#define old_ptr	(	i,
		j,
		k		)	old_ptr[i+(j+(k*ny))*nx]

Definition at line 4944 of file util_sparx.cpp.

Referenced by EMAN::Util::decimate().

#define outp	(	i,
		j,
		k		)	outp[(i+new_st_x)+((j+new_st_y)+((k+new_st_z)*new_ny))*new_nx]

Definition at line 5050 of file util_sparx.cpp.

#define outp	(	i,
		j,
		k		)	outp[i+(j+(k*new_ny))*new_nx]

Definition at line 5050 of file util_sparx.cpp.

Referenced by EMAN::Util::pad(), and EMAN::Util::window().

#define phi

(

i

)

phi [i-1]

Definition at line 6801 of file util_sparx.cpp.

Referenced by EMAN::file_store::add_image(), EMAN::OrientationGenerator::add_orientation(), ali3d_d(), EMAN::Refine3DAligner::align(), EMAN::PawelProjector::backproject3d(), EMAN::ChaoProjector::backproject3d(), EMAN::Util::even_angles(), fcalc(), fgcalc(), EMAN::RandomOrientationGenerator::gen_orientations(), EMAN::file_store::get_image(), EMAN::Transform3D::get_rotation(), EMAN::Transform::get_rotation(), EMAN::Util::hsortd(), LBD_Cart(), main(), EMAN::Util::multiref_polar_ali_2d_local(), EMAN::Util::multiref_polar_ali_2d_local_psi(), EMAN::TestImageSinewave::process_inplace(), EMAN::ChaoProjector::project3d(), EMAN::FourierGriddingProjector::project3d(), recons3d_4nn(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), refalifn3d(), EMAN::Transform3D::set_rotation(), EMAN::Transform::set_rotation(), EMAN::ChaoProjector::setdm(), slaed4_(), trans_(), unified(), EMAN::Util::vrdg(), EMAN::RT3DSphereAligner::xform_align_nbest(), and EMAN::RT3DGridAligner::xform_align_nbest().

#define PI2   QUADPI*2

Definition at line 4321 of file util_sparx.cpp.

#define PI2   2*QUADPI

Definition at line 4321 of file util_sparx.cpp.

Referenced by EMAN::Util::cml_weights(), EMAN::Util::ener(), EMAN::Util::ener_tot(), EMAN::Util::sub_fav(), and EMAN::Util::update_fav().

#define PROJ	(	i,
		j		)	PROJptr [i-1+((j-1)*NNNN)]

Definition at line 5610 of file util_sparx.cpp.

#define PROJ	(	i,
		j		)	PROJptr [i-1+((j-1)*NNNN)]

Definition at line 5610 of file util_sparx.cpp.

Referenced by EMAN::Util::WTF(), and EMAN::Util::WTM().

#define q

(

i

)

q[i-1]

Definition at line 3188 of file util_sparx.cpp.

Referenced by EMAN::Util::cluster_pairwise(), EMAN::Quaternion::create_inverse(), EMAN::Util::Crosrng_e(), crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), crosrng_ms(), EMAN::Util::Crosrng_msg(), EMAN::Util::Crosrng_msg_s(), EMAN::Util::Crosrng_msg_vec(), EMAN::Util::Crosrng_ns(), EMAN::Util::Crosrng_sm_psi(), dcstep_(), GCVmin_Tik(), EMAN::EMData::get_pixel_conv(), EMAN::EMData::get_pixel_filtered(), EMAN::Util::getBaldwinGridWeights(), inside_(), EMAN::Quaternion::interpolate(), EMAN::operator*(), EMAN::operator+(), EMAN::operator-(), EMAN::operator/(), EMAN::Util::pw_extract(), EMAN::Quaternion::Quaternion(), recons3d_CGLS_mpi_Cart(), refalin3d_perturb(), EMAN::EMData::rot_scale_conv(), EMAN::Quaternion::to_angle(), EMAN::Quaternion::to_axis(), and trfind_().

#define quadpi   3.141592653589793238462643383279502884197

Definition at line 6791 of file util_sparx.cpp.

Referenced by apmq(), and aprq2d().

#define QUADPI   3.141592653589793238462643383279502884197

Definition at line 5414 of file util_sparx.cpp.

#define QUADPI   3.141592653589793238462643383279502884197

Definition at line 5414 of file util_sparx.cpp.

#define QUADPI   3.141592653589793238462643383279502884197

Definition at line 5414 of file util_sparx.cpp.

#define rad_deg   180.0/QUADPI

Definition at line 4323 of file util_sparx.cpp.

Referenced by EMAN::Util::cml_line_in3d(), EMAN::Util::cml_line_insino(), and EMAN::Util::cml_line_insino_all().

#define rad_to_deg   180/quadpi

Definition at line 6794 of file util_sparx.cpp.

#define rad_to_dgr   180/quadpi

Definition at line 6795 of file util_sparx.cpp.

#define RI	(	i,
		j		)	RI [(i-1) + ((j-1)*3)]

Definition at line 5612 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define sign	(	x,
		y		)	(((((y)>0)?(1):(-1))*(y!=0))*(x))

Definition at line 6788 of file util_sparx.cpp.

Referenced by EMAN::Processor::EMFourierFilterFunc(), EMAN::nnSSNR_ctfReconstructor::setup(), and EMAN::nn4_ctfReconstructor::setup().

#define SS	(	I,
		J		)	SS [I-1 + (J-1)*6]

Definition at line 5611 of file util_sparx.cpp.

#define SS	(	I,
		J		)	SS [I-1 + (J-1)*6]

Definition at line 5611 of file util_sparx.cpp.

#define SS

(

I

)

SS [I-1]

Definition at line 5611 of file util_sparx.cpp.

Referenced by EMAN::Util::CANG(), EMAN::Util::WTF(), and EMAN::Util::WTM().

#define t

(

i

)

t[i-1]

Definition at line 3187 of file util_sparx.cpp.

Referenced by EMAN::OrientationGenerator::add_orientation(), EMAN::Util::ali2d_ccf_list(), EMAN::RT3DSphereAligner::align(), EMAN::RT3DGridAligner::align(), EMAN::Refine3DAligner::align(), EMAN::RefineAligner::align(), EMAN::RTFSlowExhaustiveAligner::align(), EMAN::RTFExhaustiveAligner::align(), EMAN::RotateFlipAligner::align(), EMAN::RotateTranslateFlipAligner::align(), EMAN::RotateTranslateAligner::align(), EMAN::TranslationalAligner::align(), bmv_(), EMAN::Util::BPCQ(), EMAN::Symmetry3D::cache_au_planes(), EMAN::EMData::calc_max_location(), EMAN::EMData::calc_min_location(), EMAN::EMData::calc_mutual_correlation(), EMAN::EMData::common_lines_real(), EMAN::HighpassButterworthProcessor::create_radial_func(), crlist_(), EMAN::Util::Crosrng_e(), crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), crosrng_ms(), EMAN::Util::Crosrng_msg(), EMAN::Util::Crosrng_msg_m(), EMAN::Util::Crosrng_msg_vec(), EMAN::EMData::cut_slice(), EMAN::EMData::do_radon(), EMAN::EMData::dot_rotate_translate(), dpofa_(), EMAN::Util::fftc_d(), fftc_d(), EMAN::Util::fftc_q(), fftc_q(), EMAN::Util::fftr_d(), fftr_d(), EMAN::Util::fftr_q(), fftr_q(), formk_(), EMAN::RandomOrientationGenerator::gen_orientations(), EMAN::TetrahedralSym::get_asym_unit_points(), EMAN::PlatonicSym::get_asym_unit_points(), EMAN::EMData::get_attr(), EMAN::EMData::get_pixel_filtered(), EMAN::Transform3D::get_sym_type(), EMAN::Util::get_time_label(), EMAN::Transform::icos_5_to_2(), intrsc_(), EMAN::Transform::inverse(), EMAN::Vec2< Type >::length(), EMAN::Vec3< int >::length(), main(), EMAN::Util::multiref_polar_ali_2d_local(), EMAN::Util::multiref_polar_ali_2d_local_psi(), EMAN::Transform::negate(), EMAN::FloatPoint::operator vector< float >(), EMAN::FloatSize::operator vector< float >(), EMAN::Util::point_is_in_triangle_2d(), EMAN::PawelProjector::prepcubes(), EMAN::Randnum::print_generator_type(), EMAN::ScaleTransformProcessor::process(), EMAN::TransformProcessor::process(), EMAN::TomoTiltEdgeMaskProcessor::process_inplace(), EMAN::TestTomoImage::process_inplace(), EMAN::Rotate180Processor::process_inplace(), EMAN::ScaleTransformProcessor::process_inplace(), EMAN::TransformProcessor::process_inplace(), EMAN::TestImageEllipse::process_inplace(), EMAN::TestImageHollowEllipse::process_inplace(), EMAN::IterBinMaskProcessor::process_inplace(), EMAN::AutoMask3DProcessor::process_inplace(), EMAN::SymSearchProcessor::process_inplace(), EMAN::ACFCenterProcessor::process_inplace(), EMAN::PhaseToMassCenterProcessor::process_inplace(), EMAN::ToMassCenterProcessor::process_inplace(), EMAN::FlipProcessor::process_inplace(), EMAN::NormalizeToLeastSquareProcessor::process_inplace(), EMAN::CutoffBlockProcessor::process_inplace(), EMAN::ImageProcessor::process_inplace(), EMAN::BoxMedianProcessor::process_pixel(), EMAN::StandardProjector::project3d(), refalifn(), refalifn3d(), EMAN::EMData::render_amp24(), EMAN::EMData::render_ap24(), EMAN::EMData::rot_scale_conv(), EMAN::EMData::rot_scale_conv7(), EMAN::EMData::rot_scale_conv_new(), EMAN::EMData::rot_scale_conv_new_background(), EMAN::EMData::rot_scale_trans(), EMAN::EMData::rot_scale_trans_background(), EMAN::EMData::rotate(), EMAN::Util::rotate_phase_origin(), EMAN::EMData::rotate_translate(), EMAN::Matrix4::rotation(), EMAN::EMData::scale(), EMAN::EMData::set_attr_python(), setulb_(), slaed2_(), slaed8_(), slamch_(), slasq2_(), slasq3_(), slasv2_(), sormlq_(), sormqr_(), subsm_(), EMAN::MarchingCubes::surface_face_z(), test_shared_pointer(), EMAN::Transform::tet_3_to_2(), EMAN::Gatan::to_em_datatype(), EMAN::EMData::translate(), EMAN::Transform::transpose(), trplot_(), EMAN::EMData::unwrap(), varmx(), vrplot_(), EMAN::SpiderIO::write_single_header(), EMAN::RT3DSphereAligner::xform_align_nbest(), and EMAN::RT3DGridAligner::xform_align_nbest().

#define t7

(

i

)

t7[i-1]

Definition at line 3190 of file util_sparx.cpp.

Referenced by EMAN::Util::Crosrng_e(), crosrng_e(), EMAN::Util::Crosrng_ew(), EMAN::Util::Crosrng_ms(), crosrng_ms(), EMAN::Util::Crosrng_ns(), and EMAN::Util::Crosrng_sm_psi().

#define tab1

(

i

)

tab1[i-1]

Definition at line 2639 of file util_sparx.cpp.

Referenced by EMAN::Util::fftc_d(), fftc_d(), EMAN::Util::fftc_q(), fftc_q(), EMAN::Util::fftr_d(), fftr_d(), EMAN::Util::fftr_q(), and fftr_q().

#define theta

(

i

)

theta [i-1]

Definition at line 6800 of file util_sparx.cpp.

Referenced by ali3d_d(), EMAN::PawelProjector::backproject3d(), EMAN::ChaoProjector::backproject3d(), dcstep_(), EMAN::Util::even_angles(), fcalc(), fgcalc(), EMAN::file_store::get_image(), EMAN::Util::hsortd(), LBD_Cart(), main(), mainlb_(), EMAN::Util::multiref_polar_ali_2d_local(), EMAN::Util::multiref_polar_ali_2d_local_psi(), EMAN::ChaoProjector::project3d(), EMAN::FourierGriddingProjector::project3d(), recons3d_4nn(), recons3d_CGLS_mpi_Cart(), recons3d_HyBR_mpi_Cart(), recons3d_sirt_mpi(), recons3d_sirt_mpi_Cart(), EMAN::Transform::set_rotation(), EMAN::ChaoProjector::setdm(), trans_(), unified(), and EMAN::Util::vrdg().

#define thetast

(

i

)

thetast [i-1]

Definition at line 6805 of file util_sparx.cpp.

#define TRUE   1

Definition at line 6796 of file util_sparx.cpp.

#define TRUE_   (1)

Definition at line 7500 of file util_sparx.cpp.

#define ts

(

i

)

ts [i-1]

Definition at line 6804 of file util_sparx.cpp.

#define VP

(

i

)

VP [i-1]

Definition at line 5615 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define VV

(

i

)

VV [i-1]

Definition at line 5616 of file util_sparx.cpp.

Referenced by EMAN::Util::WTM().

#define W	(	i,
		j		)	Wptr [i-1+((j-1)*Wnx)]

Definition at line 5609 of file util_sparx.cpp.

#define W	(	i,
		j		)	Wptr [i-1+((j-1)*Wnx)]

Definition at line 5609 of file util_sparx.cpp.

Referenced by EMAN::Util::getBaldwinGridWeights(), EMAN::Util::WTF(), and EMAN::Util::WTM().

#define weight

(

i

)

weight [i-1]

Definition at line 6802 of file util_sparx.cpp.

Referenced by ali3d_d(), EMAN::FRCCmp::cmp(), EMAN::FourierReconstructor::determine_slice_agreement(), EMAN::FourierReconstructor::do_insert_slice_work(), EMAN::BaldwinWoolfordReconstructor::finish(), EMAN::BaldwinWoolfordReconstructor::insert_pixel(), EMAN::BackProjectionReconstructor::insert_slice(), EMAN::WienerFourierReconstructor::insert_slice(), and EMAN::Util::vrdg().

#define xcmplx	(	i,
		j		)	xcmplx [(j-1)*2 + i-1]

Definition at line 2640 of file util_sparx.cpp.

Referenced by EMAN::Util::fftr_d(), fftr_d(), EMAN::Util::fftr_q(), and fftr_q().

#define xim	(	i,
		j		)	xim[(j-1)*nsam + i-1]

Definition at line 2161 of file util_sparx.cpp.

Referenced by EMAN::Util::bilinear(), EMAN::Util::Polar2D(), and EMAN::Util::Polar2Dm().

---

Function Documentation

int addnod_	(	int *	nst,
		int *	k,
		double *	x,
		double *	y,
		double *	z__,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	lnew,
		int *	ier
	)

Definition at line 7954 of file util_sparx.cpp.

References abs, bdyadd_(), covsph_(), intadd_(), lstptr_(), swap_(), swptst_(), and trfind_().

Referenced by trmesh_(), and EMAN::Util::trmsh3_().

{
    /* Initialized data */

    static double tol = 0.;

    /* System generated locals */
    int i__1;

    /* Local variables */
    static int l;
    static double p[3], b1, b2, b3;
    static int i1, i2, i3, kk, lp, in1, io1, io2, km1, lpf, ist, lpo1;
    extern /* Subroutine */ int swap_(int *, int *, int *,
            int *, int *, int *, int *, int *);
    static int lpo1s;
    extern /* Subroutine */ int bdyadd_(int *, int *, int *,
            int *, int *, int *, int *), intadd_(int *,
            int *, int *, int *, int *, int *, int *,
            int *), trfind_(int *, double *, int *,
            double *, double *, double *, int *, int *,
            int *, double *, double *, double *, int *,
            int *, int *), covsph_(int *, int *, int *,
            int *, int *, int *);
    extern int lstptr_(int *, int *, int *, int *);
    extern long int swptst_(int *, int *, int *, int *,
            double *, double *, double *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   01/08/03 */

/*   This subroutine adds node K to a triangulation of the */
/* convex hull of nodes 1,...,K-1, producing a triangulation */
/* of the convex hull of nodes 1,...,K. */

/*   The algorithm consists of the following steps:  node K */
/* is located relative to the triangulation (TRFIND), its */
/* index is added to the data structure (INTADD or BDYADD), */
/* and a sequence of swaps (SWPTST and SWAP) are applied to */
/* the arcs opposite K so that all arcs incident on node K */
/* and opposite node K are locally optimal (satisfy the cir- */
/* cumcircle test).  Thus, if a Delaunay triangulation is */
/* input, a Delaunay triangulation will result. */


/* On input: */

/*       NST = Index of a node at which TRFIND begins its */
/*             search.  Search time depends on the proximity */
/*             of this node to K.  If NST < 1, the search is */
/*             begun at node K-1. */

/*       K = Nodal index (index for X, Y, Z, and LEND) of the */
/*           new node to be added.  K .GE. 4. */

/*       X,Y,Z = Arrays of length .GE. K containing Car- */
/*               tesian coordinates of the nodes. */
/*               (X(I),Y(I),Z(I)) defines node I for */
/*               I = 1,...,K. */

/* The above parameters are not altered by this routine. */

/*       LIST,LPTR,LEND,LNEW = Data structure associated with */
/*                             the triangulation of nodes 1 */
/*                             to K-1.  The array lengths are */
/*                             assumed to be large enough to */
/*                             add node K.  Refer to Subrou- */
/*                             tine TRMESH. */

/* On output: */

/*       LIST,LPTR,LEND,LNEW = Data structure updated with */
/*                             the addition of node K as the */
/*                             last entry unless IER .NE. 0 */
/*                             and IER .NE. -3, in which case */
/*                             the arrays are not altered. */

/*       IER = Error indicator: */
/*             IER =  0 if no errors were encountered. */
/*             IER = -1 if K is outside its valid range */
/*                      on input. */
/*             IER = -2 if all nodes (including K) are col- */
/*                      linear (lie on a common geodesic). */
/*             IER =  L if nodes L and K coincide for some */
/*                      L < K.  Refer to TOL below. */

/* Modules required by ADDNOD:  BDYADD, COVSPH, INSERT, */
/*                                INTADD, JRAND, LSTPTR, */
/*                                STORE, SWAP, SWPTST, */
/*                                TRFIND */

/* Intrinsic function called by ADDNOD:  ABS */

/* *********************************************************** */


/* Local parameters: */

/* B1,B2,B3 = Unnormalized barycentric coordinates returned */
/*              by TRFIND. */
/* I1,I2,I3 = Vertex indexes of a triangle containing K */
/* IN1 =      Vertex opposite K:  first neighbor of IO2 */
/*              that precedes IO1.  IN1,IO1,IO2 are in */
/*              counterclockwise order. */
/* IO1,IO2 =  Adjacent neighbors of K defining an arc to */
/*              be tested for a swap */
/* IST =      Index of node at which TRFIND begins its search */
/* KK =       Local copy of K */
/* KM1 =      K-1 */
/* L =        Vertex index (I1, I2, or I3) returned in IER */
/*              if node K coincides with a vertex */
/* LP =       LIST pointer */
/* LPF =      LIST pointer to the first neighbor of K */
/* LPO1 =     LIST pointer to IO1 */
/* LPO1S =    Saved value of LPO1 */
/* P =        Cartesian coordinates of node K */
/* TOL =      Tolerance defining coincident nodes:  bound on */
/*              the deviation from 1 of the cosine of the */
/*              angle between the nodes. */
/*              Note that |1-cos(A)| is approximately A*A/2. */

    /* Parameter adjustments */
    --lend;
    --z__;
    --y;
    --x;
    --list;
    --lptr;

    /* Function Body */

    kk = *k;
    if (kk < 4) {
        goto L3;
    }

/* Initialization: */
    km1 = kk - 1;
    ist = *nst;
    if (ist < 1) {
        ist = km1;
    }
    p[0] = x[kk];
    p[1] = y[kk];
    p[2] = z__[kk];

/* Find a triangle (I1,I2,I3) containing K or the rightmost */
/*   (I1) and leftmost (I2) visible boundary nodes as viewed */
/*   from node K. */
    trfind_(&ist, p, &km1, &x[1], &y[1], &z__[1], &list[1], &lptr[1], &lend[1]
            , &b1, &b2, &b3, &i1, &i2, &i3);

/*   Test for collinear or (nearly) duplicate nodes. */

    if (i1 == 0) {
        goto L4;
    }
    l = i1;
    if (p[0] * x[l] + p[1] * y[l] + p[2] * z__[l] >= 1. - tol) {
        goto L5;
    }
    l = i2;
    if (p[0] * x[l] + p[1] * y[l] + p[2] * z__[l] >= 1. - tol) {
        goto L5;
    }
    if (i3 != 0) {
        l = i3;
        if (p[0] * x[l] + p[1] * y[l] + p[2] * z__[l] >= 1. - tol) {
            goto L5;
        }
        intadd_(&kk, &i1, &i2, &i3, &list[1], &lptr[1], &lend[1], lnew);
    } else {
        if (i1 != i2) {
            bdyadd_(&kk, &i1, &i2, &list[1], &lptr[1], &lend[1], lnew);
        } else {
            covsph_(&kk, &i1, &list[1], &lptr[1], &lend[1], lnew);
        }
    }
    *ier = 0;

/* Initialize variables for optimization of the */
/*   triangulation. */
    lp = lend[kk];
    lpf = lptr[lp];
    io2 = list[lpf];
    lpo1 = lptr[lpf];
    io1 = (i__1 = list[lpo1], abs(i__1));

/* Begin loop:  find the node opposite K. */

L1:
    lp = lstptr_(&lend[io1], &io2, &list[1], &lptr[1]);
    if (list[lp] < 0) {
        goto L2;
    }
    lp = lptr[lp];
    in1 = (i__1 = list[lp], abs(i__1));

/* Swap test:  if a swap occurs, two new arcs are */
/*             opposite K and must be tested. */

    lpo1s = lpo1;
    if (! swptst_(&in1, &kk, &io1, &io2, &x[1], &y[1], &z__[1])) {
        goto L2;
    }
    swap_(&in1, &kk, &io1, &io2, &list[1], &lptr[1], &lend[1], &lpo1);
    if (lpo1 == 0) {

/*   A swap is not possible because KK and IN1 are already */
/*     adjacent.  This error in SWPTST only occurs in the */
/*     neutral case and when there are nearly duplicate */
/*     nodes. */

        lpo1 = lpo1s;
        goto L2;
    }
    io1 = in1;
    goto L1;

/* No swap occurred.  Test for termination and reset */
/*   IO2 and IO1. */

L2:
    if (lpo1 == lpf || list[lpo1] < 0) {
        return 0;
    }
    io2 = io1;
    lpo1 = lptr[lpo1];
    io1 = (i__1 = list[lpo1], abs(i__1));
    goto L1;

/* KK < 4. */

L3:
    *ier = -1;
    return 0;

/* All nodes are collinear. */

L4:
    *ier = -2;
    return 0;

/* Nodes L and K coincide. */

L5:
    *ier = l;
    return 0;
} /* addnod_ */

double angle_	(	double *	v1,
		double *	v2,
		double *	v3
	)

Definition at line 8213 of file util_sparx.cpp.

References left_(), and sqrt().

Referenced by areav_new__().

{
    /* System generated locals */
    double ret_val;

    /* Builtin functions */
    //double sqrt(double), acos(double);

    /* Local variables */
    static double a;
    static int i__;
    static double ca, s21, s23, u21[3], u23[3];
    extern long int left_(double *, double *, double *, double
            *, double *, double *, double *, double *,
            double *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   06/03/03 */

/*   Given a sequence of three nodes (V1,V2,V3) on the sur- */
/* face of the unit sphere, this function returns the */
/* interior angle at V2 -- the dihedral angle between the */
/* plane defined by V2 and V3 (and the origin) and the plane */
/* defined by V2 and V1 or, equivalently, the angle between */
/* the normals V2 X V3 and V2 X V1.  Note that the angle is */
/* in the range 0 to Pi if V3 Left V1->V2, Pi to 2*Pi other- */
/* wise.  The surface area of a spherical polygon with CCW- */
/* ordered vertices V1, V2, ..., Vm is Asum - (m-2)*Pi, where */
/* Asum is the sum of the m interior angles computed from the */
/* sequences (Vm,V1,V2), (V1,V2,V3), (V2,V3,V4), ..., */
/* (Vm-1,Vm,V1). */


/* On input: */

/*       V1,V2,V3 = Arrays of length 3 containing the Carte- */
/*                  sian coordinates of unit vectors.  These */
/*                  vectors, if nonzero, are implicitly */
/*                  scaled to have length 1. */

/* Input parameters are not altered by this function. */

/* On output: */

/*       ANGLE = Angle defined above, or 0 if V2 X V1 = 0 or */
/*               V2 X V3 = 0. */

/* Module required by ANGLE:  LEFT */

/* Intrinsic functions called by ANGLE:  ACOS, SQRT */

/* *********************************************************** */


/* Local parameters: */

/* A =       Interior angle at V2 */
/* CA =      cos(A) */
/* I =       DO-loop index and index for U21 and U23 */
/* S21,S23 = Sum of squared components of U21 and U23 */
/* U21,U23 = Unit normal vectors to the planes defined by */
/*             pairs of triangle vertices */


/* Compute cross products U21 = V2 X V1 and U23 = V2 X V3. */

    /* Parameter adjustments */
    --v3;
    --v2;
    --v1;

    /* Function Body */
    u21[0] = v2[2] * v1[3] - v2[3] * v1[2];
    u21[1] = v2[3] * v1[1] - v2[1] * v1[3];
    u21[2] = v2[1] * v1[2] - v2[2] * v1[1];

    u23[0] = v2[2] * v3[3] - v2[3] * v3[2];
    u23[1] = v2[3] * v3[1] - v2[1] * v3[3];
    u23[2] = v2[1] * v3[2] - v2[2] * v3[1];

/* Normalize U21 and U23 to unit vectors. */

    s21 = 0.;
    s23 = 0.;
    for (i__ = 1; i__ <= 3; ++i__) {
        s21 += u21[i__ - 1] * u21[i__ - 1];
        s23 += u23[i__ - 1] * u23[i__ - 1];
/* L1: */
    }

/* Test for a degenerate triangle associated with collinear */
/*   vertices. */

    if (s21 == 0. || s23 == 0.) {
        ret_val = 0.;
        return ret_val;
    }
    s21 = sqrt(s21);
    s23 = sqrt(s23);
    for (i__ = 1; i__ <= 3; ++i__) {
        u21[i__ - 1] /= s21;
        u23[i__ - 1] /= s23;
/* L2: */
    }

/* Compute the angle A between normals: */

/*   CA = cos(A) = <U21,U23> */

    ca = u21[0] * u23[0] + u21[1] * u23[1] + u21[2] * u23[2];
    if (ca < -1.) {
        ca = -1.;
    }
    if (ca > 1.) {
        ca = 1.;
    }
    a = acos(ca);

/* Adjust A to the interior angle:  A > Pi iff */
/*   V3 Right V1->V2. */

    if (! left_(&v1[1], &v1[2], &v1[3], &v2[1], &v2[2], &v2[3], &v3[1], &v3[2]
            , &v3[3])) {
        a = acos(-1.) * 2. - a;
    }
    ret_val = a;
    return ret_val;
} /* angle_ */

double areas_	(	double *	v1,
		double *	v2,
		double *	v3
	)

Definition at line 8349 of file util_sparx.cpp.

References sqrt().

Referenced by EMAN::Util::areav_().

{
    /* System generated locals */
    double ret_val;

    /* Builtin functions */
    //double sqrt(double), acos(double);

    /* Local variables */
    static int i__;
    static double a1, a2, a3, s12, s31, s23, u12[3], u23[3], u31[3], ca1,
            ca2, ca3;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   06/22/98 */

/*   This function returns the area of a spherical triangle */
/* on the unit sphere. */


/* On input: */

/*       V1,V2,V3 = Arrays of length 3 containing the Carte- */
/*                  sian coordinates of unit vectors (the */
/*                  three triangle vertices in any order). */
/*                  These vectors, if nonzero, are implicitly */
/*                  scaled to have length 1. */

/* Input parameters are not altered by this function. */

/* On output: */

/*       AREAS = Area of the spherical triangle defined by */
/*               V1, V2, and V3 in the range 0 to 2*PI (the */
/*               area of a hemisphere).  AREAS = 0 (or 2*PI) */
/*               if and only if V1, V2, and V3 lie in (or */
/*               close to) a plane containing the origin. */

/* Modules required by AREAS:  None */

/* Intrinsic functions called by AREAS:  ACOS, SQRT */

/* *********************************************************** */


/* Local parameters: */

/* A1,A2,A3 =    Interior angles of the spherical triangle */
/* CA1,CA2,CA3 = cos(A1), cos(A2), and cos(A3), respectively */
/* I =           DO-loop index and index for Uij */
/* S12,S23,S31 = Sum of squared components of U12, U23, U31 */
/* U12,U23,U31 = Unit normal vectors to the planes defined by */
/*                 pairs of triangle vertices */


/* Compute cross products Uij = Vi X Vj. */

    /* Parameter adjustments */
    --v3;
    --v2;
    --v1;

    /* Function Body */
    u12[0] = v1[2] * v2[3] - v1[3] * v2[2];
    u12[1] = v1[3] * v2[1] - v1[1] * v2[3];
    u12[2] = v1[1] * v2[2] - v1[2] * v2[1];

    u23[0] = v2[2] * v3[3] - v2[3] * v3[2];
    u23[1] = v2[3] * v3[1] - v2[1] * v3[3];
    u23[2] = v2[1] * v3[2] - v2[2] * v3[1];

    u31[0] = v3[2] * v1[3] - v3[3] * v1[2];
    u31[1] = v3[3] * v1[1] - v3[1] * v1[3];
    u31[2] = v3[1] * v1[2] - v3[2] * v1[1];

/* Normalize Uij to unit vectors. */

    s12 = 0.;
    s23 = 0.;
    s31 = 0.;
    for (i__ = 1; i__ <= 3; ++i__) {
        s12 += u12[i__ - 1] * u12[i__ - 1];
        s23 += u23[i__ - 1] * u23[i__ - 1];
        s31 += u31[i__ - 1] * u31[i__ - 1];
/* L2: */
    }

/* Test for a degenerate triangle associated with collinear */
/*   vertices. */

    if (s12 == 0. || s23 == 0. || s31 == 0.) {
        ret_val = 0.;
        return ret_val;
    }
    s12 = sqrt(s12);
    s23 = sqrt(s23);
    s31 = sqrt(s31);
    for (i__ = 1; i__ <= 3; ++i__) {
        u12[i__ - 1] /= s12;
        u23[i__ - 1] /= s23;
        u31[i__ - 1] /= s31;
/* L3: */
    }

/* Compute interior angles Ai as the dihedral angles between */
/*   planes: */
/*           CA1 = cos(A1) = -<U12,U31> */
/*           CA2 = cos(A2) = -<U23,U12> */
/*           CA3 = cos(A3) = -<U31,U23> */

    ca1 = -u12[0] * u31[0] - u12[1] * u31[1] - u12[2] * u31[2];
    ca2 = -u23[0] * u12[0] - u23[1] * u12[1] - u23[2] * u12[2];
    ca3 = -u31[0] * u23[0] - u31[1] * u23[1] - u31[2] * u23[2];
    if (ca1 < -1.) {
        ca1 = -1.;
    }
    if (ca1 > 1.) {
        ca1 = 1.;
    }
    if (ca2 < -1.) {
        ca2 = -1.;
    }
    if (ca2 > 1.) {
        ca2 = 1.;
    }
    if (ca3 < -1.) {
        ca3 = -1.;
    }
    if (ca3 > 1.) {
        ca3 = 1.;
    }
    a1 = acos(ca1);
    a2 = acos(ca2);
    a3 = acos(ca3);

/* Compute AREAS = A1 + A2 + A3 - PI. */

    ret_val = a1 + a2 + a3 - acos(-1.);
    if (ret_val < 0.) {
        ret_val = 0.;
    }
    return ret_val;
} /* areas_ */

double areav_new__	(	int *	k,
		int *	n,
		double *	x,
		double *	y,
		double *	z__,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	ier
	)

Definition at line 8704 of file util_sparx.cpp.

References angle_(), circum_(), and ierr.

{
    /* System generated locals */
    double ret_val = 0;

    /* Builtin functions */
    //double acos(double);

    /* Local variables */
    static int m;
    static double c1[3], c2[3], c3[3];
    static int n1, n2, n3;
    static double v1[3], v2[3], v3[3];
    static int lp;
    static double c1s[3], c2s[3];
    static int lpl, ierr;
    static double asum;
    extern double angle_(double *, double *, double *);
    static float areav;
    extern /* Subroutine */ int circum_(double *, double *,
            double *, double *, int *);


/* *********************************************************** */

/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   06/03/03 */

/*   Given a Delaunay triangulation and the index K of an */
/* interior node, this subroutine returns the (surface) area */
/* of the Voronoi region associated with node K.  The Voronoi */
/* region is the polygon whose vertices are the circumcenters */
/* of the triangles that contain node K, where a triangle */
/* circumcenter is the point (unit vector) lying at the same */
/* angular distance from the three vertices and contained in */
/* the same hemisphere as the vertices.  The Voronoi region */
/* area is computed as Asum-(m-2)*Pi, where m is the number */
/* of Voronoi vertices (neighbors of K) and Asum is the sum */
/* of interior angles at the vertices. */


/* On input: */

/*       K = Nodal index in the range 1 to N. */

/*       N = Number of nodes in the triangulation.  N > 3. */

/*       X,Y,Z = Arrays of length N containing the Cartesian */
/*               coordinates of the nodes (unit vectors). */

/*       LIST,LPTR,LEND = Data structure defining the trian- */
/*                        gulation.  Refer to Subroutine */
/*                        TRMESH. */

/* Input parameters are not altered by this function. */

/* On output: */

/*       AREAV = Area of Voronoi region K unless IER > 0, */
/*               in which case AREAV = 0. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if K or N is outside its valid range */
/*                     on input. */
/*             IER = 2 if K indexes a boundary node. */
/*             IER = 3 if an error flag is returned by CIRCUM */
/*                     (null triangle). */

/* Modules required by AREAV:  ANGLE, CIRCUM */

/* Intrinsic functions called by AREAV:  ACOS, DBLE */

/* *********************************************************** */


/* Test for invalid input. */

    /* Parameter adjustments */
    --lend;
    --z__;
    --y;
    --x;
    --list;
    --lptr;

    /* Function Body */
    if (*k < 1 || *k > *n || *n <= 3) {
        goto L11;
    }

/* Initialization:  Set N3 to the last neighbor of N1 = K. */
/*   The number of neighbors and the sum of interior angles */
/*   are accumulated in M and ASUM, respectively. */

    n1 = *k;
    v1[0] = x[n1];
    v1[1] = y[n1];
    v1[2] = z__[n1];
    lpl = lend[n1];
    n3 = list[lpl];
    if (n3 < 0) {
        goto L12;
    }
    lp = lpl;
    m = 0;
    asum = 0.;

/* Loop on triangles (N1,N2,N3) containing N1 = K. */

L1:
    ++m;
    n2 = n3;
    lp = lptr[lp];
    n3 = list[lp];
    v2[0] = x[n2];
    v2[1] = y[n2];
    v2[2] = z__[n2];
    v3[0] = x[n3];
    v3[1] = y[n3];
    v3[2] = z__[n3];
    if (m == 1) {

/* First triangle:  compute the circumcenter C2 and save a */
/*   copy in C1S. */

        circum_(v1, v2, v3, c2, &ierr);
        if (ierr != 0) {
            goto L13;
        }
        c1s[0] = c2[0];
        c1s[1] = c2[1];
        c1s[2] = c2[2];
    } else if (m == 2) {

/* Second triangle:  compute the circumcenter C3 and save a */
/*   copy in C2S. */

        circum_(v1, v2, v3, c3, &ierr);
        if (ierr != 0) {
            goto L13;
        }
        c2s[0] = c3[0];
        c2s[1] = c3[1];
        c2s[2] = c3[2];
    } else {

/* Set C1 to C2, set C2 to C3, compute the new circumcenter */
/*   C3, and compute the interior angle at C2 from the */
/*   sequence of vertices (C1,C2,C3). */

        c1[0] = c2[0];
        c1[1] = c2[1];
        c1[2] = c2[2];
        c2[0] = c3[0];
        c2[1] = c3[1];
        c2[2] = c3[2];
        circum_(v1, v2, v3, c3, &ierr);
        if (ierr != 0) {
            goto L13;
        }
        asum += angle_(c1, c2, c3);
    }

/* Bottom on loop on neighbors of K. */

    if (lp != lpl) {
        goto L1;
    }

/* C3 is the last vertex.  Compute its interior angle from */
/*   the sequence (C2,C3,C1S). */

    asum += angle_(c2, c3, c1s);

/* Compute the interior angle at C1S from */
/*   the sequence (C3,C1S,C2S). */

    asum += angle_(c3, c1s, c2s);

/* No error encountered. */

    *ier = 0;
    ret_val = asum - (double) (m - 2) * acos(-1.);
    return ret_val;

/* Invalid input. */

L11:
    *ier = 1;
    areav = 0.f;
    return ret_val;

/* K indexes a boundary node. */

L12:
    *ier = 2;
    areav = 0.f;
    return ret_val;

/* Error in CIRCUM. */

L13:
    *ier = 3;
    areav = 0.f;
    return ret_val;
} /* areav_new__ */

int bdyadd_	(	int *	kk,
		int *	i1,
		int *	i2,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	lnew
	)

Definition at line 8917 of file util_sparx.cpp.

References insert_().

Referenced by addnod_().

{
    static int k, n1, n2, lp, lsav, nsav, next;
    extern /* Subroutine */ int insert_(int *, int *, int *,
            int *, int *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/11/96 */

/*   This subroutine adds a boundary node to a triangulation */
/* of a set of KK-1 points on the unit sphere.  The data */
/* structure is updated with the insertion of node KK, but no */
/* optimization is performed. */

/*   This routine is identical to the similarly named routine */
/* in TRIPACK. */


/* On input: */

/*       KK = Index of a node to be connected to the sequence */
/*            of all visible boundary nodes.  KK .GE. 1 and */
/*            KK must not be equal to I1 or I2. */

/*       I1 = First (rightmost as viewed from KK) boundary */
/*            node in the triangulation that is visible from */
/*            node KK (the line segment KK-I1 intersects no */
/*            arcs. */

/*       I2 = Last (leftmost) boundary node that is visible */
/*            from node KK.  I1 and I2 may be determined by */
/*            Subroutine TRFIND. */

/* The above parameters are not altered by this routine. */

/*       LIST,LPTR,LEND,LNEW = Triangulation data structure */
/*                             created by Subroutine TRMESH. */
/*                             Nodes I1 and I2 must be in- */
/*                             cluded in the triangulation. */

/* On output: */

/*       LIST,LPTR,LEND,LNEW = Data structure updated with */
/*                             the addition of node KK.  Node */
/*                             KK is connected to I1, I2, and */
/*                             all boundary nodes in between. */

/* Module required by BDYADD:  INSERT */

/* *********************************************************** */


/* Local parameters: */

/* K =     Local copy of KK */
/* LP =    LIST pointer */
/* LSAV =  LIST pointer */
/* N1,N2 = Local copies of I1 and I2, respectively */
/* NEXT =  Boundary node visible from K */
/* NSAV =  Boundary node visible from K */

    /* Parameter adjustments */
    --lend;
    --lptr;
    --list;

    /* Function Body */
    k = *kk;
    n1 = *i1;
    n2 = *i2;

/* Add K as the last neighbor of N1. */

    lp = lend[n1];
    lsav = lptr[lp];
    lptr[lp] = *lnew;
    list[*lnew] = -k;
    lptr[*lnew] = lsav;
    lend[n1] = *lnew;
    ++(*lnew);
    next = -list[lp];
    list[lp] = next;
    nsav = next;

/* Loop on the remaining boundary nodes between N1 and N2, */
/*   adding K as the first neighbor. */

L1:
    lp = lend[next];
    insert_(&k, &lp, &list[1], &lptr[1], lnew);
    if (next == n2) {
        goto L2;
    }
    next = -list[lp];
    list[lp] = next;
    goto L1;

/* Add the boundary nodes between N1 and N2 as neighbors */
/*   of node K. */

L2:
    lsav = *lnew;
    list[*lnew] = n1;
    lptr[*lnew] = *lnew + 1;
    ++(*lnew);
    next = nsav;

L3:
    if (next == n2) {
        goto L4;
    }
    list[*lnew] = next;
    lptr[*lnew] = *lnew + 1;
    ++(*lnew);
    lp = lend[next];
    next = list[lp];
    goto L3;

L4:
    list[*lnew] = -n2;
    lptr[*lnew] = lsav;
    lend[k] = *lnew;
    ++(*lnew);
    return 0;
} /* bdyadd_ */

int bnodes_	(	int *	n,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	nodes,
		int *	nb,
		int *	na,
		int *	nt
	)

Definition at line 9051 of file util_sparx.cpp.

References nn().

{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int k, n0, lp, nn, nst;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   06/26/96 */

/*   Given a triangulation of N nodes on the unit sphere */
/* created by Subroutine TRMESH, this subroutine returns an */
/* array containing the indexes (if any) of the counterclock- */
/* wise-ordered sequence of boundary nodes -- the nodes on */
/* the boundary of the convex hull of the set of nodes.  (The */
/* boundary is empty if the nodes do not lie in a single */
/* hemisphere.)  The numbers of boundary nodes, arcs, and */
/* triangles are also returned. */


/* On input: */

/*       N = Number of nodes in the triangulation.  N .GE. 3. */

/*       LIST,LPTR,LEND = Data structure defining the trian- */
/*                        gulation.  Refer to Subroutine */
/*                        TRMESH. */

/* The above parameters are not altered by this routine. */

/*       NODES = int array of length at least NB */
/*               (NB .LE. N). */

/* On output: */

/*       NODES = Ordered sequence of boundary node indexes */
/*               in the range 1 to N (in the first NB loca- */
/*               tions). */

/*       NB = Number of boundary nodes. */

/*       NA,NT = Number of arcs and triangles, respectively, */
/*               in the triangulation. */

/* Modules required by BNODES:  None */

/* *********************************************************** */


/* Local parameters: */

/* K =   NODES index */
/* LP =  LIST pointer */
/* N0 =  Boundary node to be added to NODES */
/* NN =  Local copy of N */
/* NST = First element of nodes (arbitrarily chosen to be */
/*         the one with smallest index) */

    /* Parameter adjustments */
    --lend;
    --list;
    --lptr;
    --nodes;

    /* Function Body */
    nn = *n;

/* Search for a boundary node. */

    i__1 = nn;
    for (nst = 1; nst <= i__1; ++nst) {
        lp = lend[nst];
        if (list[lp] < 0) {
            goto L2;
        }
/* L1: */
    }

/* The triangulation contains no boundary nodes. */

    *nb = 0;
    *na = (nn - 2) * 3;
    *nt = nn - (2<<1);
    return 0;

/* NST is the first boundary node encountered.  Initialize */
/*   for traversal of the boundary. */

L2:
    nodes[1] = nst;
    k = 1;
    n0 = nst;

/* Traverse the boundary in counterclockwise order. */

L3:
    lp = lend[n0];
    lp = lptr[lp];
    n0 = list[lp];
    if (n0 == nst) {
        goto L4;
    }
    ++k;
    nodes[k] = n0;
    goto L3;

/* Store the counts. */

L4:
    *nb = k;
    *nt = (*n << 1) - *nb - 2;
    *na = *nt + *n - 1;
    return 0;
} /* bnodes_ */

int circle_	(	int *	k,
		double *	xc,
		double *	yc,
		int *	ier
	)

Definition at line 9175 of file util_sparx.cpp.

{
    /* System generated locals */
    int i__1;

    /* Builtin functions */
    //double atan(double), cos(double), sin(double);

    /* Local variables */
    static double a, c__;
    static int i__;
    static double s;
    static int k2, k3;
    static double x0, y0;
    static int kk, np1;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   04/06/90 */

/*   This subroutine computes the coordinates of a sequence */
/* of N equally spaced points on the unit circle centered at */
/* (0,0).  An N-sided polygonal approximation to the circle */
/* may be plotted by connecting (XC(I),YC(I)) to (XC(I+1), */
/* YC(I+1)) for I = 1,...,N, where XC(N+1) = XC(1) and */
/* YC(N+1) = YC(1).  A reasonable value for N in this case */
/* is 2*PI*R, where R is the radius of the circle in device */
/* coordinates. */


/* On input: */

/*       K = Number of points in each quadrant, defining N as */
/*           4K.  K .GE. 1. */

/*       XC,YC = Arrays of length at least N+1 = 4K+1. */

/* K is not altered by this routine. */

/* On output: */

/*       XC,YC = Cartesian coordinates of the points on the */
/*               unit circle in the first N+1 locations. */
/*               XC(I) = cos(A*(I-1)), YC(I) = sin(A*(I-1)), */
/*               where A = 2*PI/N.  Note that XC(N+1) = XC(1) */
/*               and YC(N+1) = YC(1). */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if K < 1 on input. */

/* Modules required by CIRCLE:  None */

/* Intrinsic functions called by CIRCLE:  ATAN, COS, DBLE, */
/*                                          SIN */

/* *********************************************************** */


/* Local parameters: */

/* I =     DO-loop index and index for XC and YC */
/* KK =    Local copy of K */
/* K2 =    K*2 */
/* K3 =    K*3 */
/* NP1 =   N+1 = 4*K + 1 */
/* A =     Angular separation between adjacent points */
/* C,S =   Cos(A) and sin(A), respectively, defining a */
/*           rotation through angle A */
/* X0,Y0 = Cartesian coordinates of a point on the unit */
/*           circle in the first quadrant */

    /* Parameter adjustments */
    --yc;
    --xc;

    /* Function Body */
    kk = *k;
    k2 = kk << 1;
    k3 = kk * 3;
    np1 = (kk << 2) + 1;

/* Test for invalid input, compute A, C, and S, and */
/*   initialize (X0,Y0) to (1,0). */

    if (kk < 1) {
        goto L2;
    }
    a = atan(1.) * 2. / (double) kk;
    c__ = cos(a);
    s = sin(a);
    x0 = 1.;
    y0 = 0.;

/* Loop on points (X0,Y0) in the first quadrant, storing */
/*   the point and its reflections about the x axis, the */
/*   y axis, and the line y = -x. */

    i__1 = kk;
    for (i__ = 1; i__ <= i__1; ++i__) {
        xc[i__] = x0;
        yc[i__] = y0;
        xc[i__ + kk] = -y0;
        yc[i__ + kk] = x0;
        xc[i__ + k2] = -x0;
        yc[i__ + k2] = -y0;
        xc[i__ + k3] = y0;
        yc[i__ + k3] = -x0;

/*   Rotate (X0,Y0) counterclockwise through angle A. */

        x0 = c__ * x0 - s * y0;
        y0 = s * x0 + c__ * y0;
/* L1: */
    }

/* Store the coordinates of the first point as the last */
/*   point. */

    xc[np1] = xc[1];
    yc[np1] = yc[1];
    *ier = 0;
    return 0;

/* K < 1. */

L2:
    *ier = 1;
    return 0;
} /* circle_ */

int circum_	(	double *	v1,
		double *	v2,
		double *	v3,
		double *	c__,
		int *	ier
	)

Definition at line 9313 of file util_sparx.cpp.

References sqrt().

Referenced by EMAN::Util::areav_(), areav_new__(), and crlist_().

{
    /* Builtin functions */
    //double sqrt(double);

    /* Local variables */
    static int i__;
    static double e1[3], e2[3], cu[3], cnorm;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   10/27/02 */

/*   This subroutine returns the circumcenter of a spherical */
/* triangle on the unit sphere:  the point on the sphere sur- */
/* face that is equally distant from the three triangle */
/* vertices and lies in the same hemisphere, where distance */
/* is taken to be arc-length on the sphere surface. */


/* On input: */

/*       V1,V2,V3 = Arrays of length 3 containing the Carte- */
/*                  sian coordinates of the three triangle */
/*                  vertices (unit vectors) in CCW order. */

/* The above parameters are not altered by this routine. */

/*       C = Array of length 3. */

/* On output: */

/*       C = Cartesian coordinates of the circumcenter unless */
/*           IER > 0, in which case C is not defined.  C = */
/*           (V2-V1) X (V3-V1) normalized to a unit vector. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if V1, V2, and V3 lie on a common */
/*                     line:  (V2-V1) X (V3-V1) = 0. */
/*             (The vertices are not tested for validity.) */

/* Modules required by CIRCUM:  None */

/* Intrinsic function called by CIRCUM:  SQRT */

/* *********************************************************** */


/* Local parameters: */

/* CNORM = Norm of CU:  used to compute C */
/* CU =    Scalar multiple of C:  E1 X E2 */
/* E1,E2 = Edges of the underlying planar triangle: */
/*           V2-V1 and V3-V1, respectively */
/* I =     DO-loop index */

    /* Parameter adjustments */
    --c__;
    --v3;
    --v2;
    --v1;

    /* Function Body */
    for (i__ = 1; i__ <= 3; ++i__) {
        e1[i__ - 1] = v2[i__] - v1[i__];
        e2[i__ - 1] = v3[i__] - v1[i__];
/* L1: */
    }

/* Compute CU = E1 X E2 and CNORM**2. */

    cu[0] = e1[1] * e2[2] - e1[2] * e2[1];
    cu[1] = e1[2] * e2[0] - e1[0] * e2[2];
    cu[2] = e1[0] * e2[1] - e1[1] * e2[0];
    cnorm = cu[0] * cu[0] + cu[1] * cu[1] + cu[2] * cu[2];

/* The vertices lie on a common line if and only if CU is */
/*   the zero vector. */

    if (cnorm != 0.) {

/*   No error:  compute C. */

        cnorm = sqrt(cnorm);
        for (i__ = 1; i__ <= 3; ++i__) {
            c__[i__] = cu[i__ - 1] / cnorm;
/* L2: */
        }

/* If the vertices are nearly identical, the problem is */
/*   ill-conditioned and it is possible for the computed */
/*   value of C to be 180 degrees off:  <C,V1> near -1 */
/*   when it should be positive. */

        if (c__[1] * v1[1] + c__[2] * v1[2] + c__[3] * v1[3] < -.5) {
            c__[1] = -c__[1];
            c__[2] = -c__[2];
            c__[3] = -c__[3];
        }
        *ier = 0;
    } else {

/*   CU = 0. */

        *ier = 1;
    }
    return 0;
} /* circum_ */

int covsph_	(	int *	kk,
		int *	n0,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	lnew
	)

Definition at line 9430 of file util_sparx.cpp.

References insert_().

Referenced by addnod_().

{
    static int k, lp, nst, lsav, next;
    extern /* Subroutine */ int insert_(int *, int *, int *,
            int *, int *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/17/96 */

/*   This subroutine connects an exterior node KK to all */
/* boundary nodes of a triangulation of KK-1 points on the */
/* unit sphere, producing a triangulation that covers the */
/* sphere.  The data structure is updated with the addition */
/* of node KK, but no optimization is performed.  All boun- */
/* dary nodes must be visible from node KK. */


/* On input: */

/*       KK = Index of the node to be connected to the set of */
/*            all boundary nodes.  KK .GE. 4. */

/*       N0 = Index of a boundary node (in the range 1 to */
/*            KK-1).  N0 may be determined by Subroutine */
/*            TRFIND. */

/* The above parameters are not altered by this routine. */

/*       LIST,LPTR,LEND,LNEW = Triangulation data structure */
/*                             created by Subroutine TRMESH. */
/*                             Node N0 must be included in */
/*                             the triangulation. */

/* On output: */

/*       LIST,LPTR,LEND,LNEW = Data structure updated with */
/*                             the addition of node KK as the */
/*                             last entry.  The updated */
/*                             triangulation contains no */
/*                             boundary nodes. */

/* Module required by COVSPH:  INSERT */

/* *********************************************************** */


/* Local parameters: */

/* K =     Local copy of KK */
/* LP =    LIST pointer */
/* LSAV =  LIST pointer */
/* NEXT =  Boundary node visible from K */
/* NST =   Local copy of N0 */

    /* Parameter adjustments */
    --lend;
    --lptr;
    --list;

    /* Function Body */
    k = *kk;
    nst = *n0;

/* Traverse the boundary in clockwise order, inserting K as */
/*   the first neighbor of each boundary node, and converting */
/*   the boundary node to an interior node. */

    next = nst;
L1:
    lp = lend[next];
    insert_(&k, &lp, &list[1], &lptr[1], lnew);
    next = -list[lp];
    list[lp] = next;
    if (next != nst) {
        goto L1;
    }

/* Traverse the boundary again, adding each node to K's */
/*   adjacency list. */

    lsav = *lnew;
L2:
    lp = lend[next];
    list[*lnew] = next;
    lptr[*lnew] = *lnew + 1;
    ++(*lnew);
    next = list[lp];
    if (next != nst) {
        goto L2;
    }

    lptr[*lnew - 1] = lsav;
    lend[k] = *lnew - 1;
    return 0;
} /* covsph_ */

int crlist_	(	int *	n,
		int *	ncol,
		double *	x,
		double *	y,
		double *	z__,
		int *	list,
		int *	lend,
		int *	lptr,
		int *	lnew,
		int *	ltri,
		int *	listc,
		int *	nb,
		double *	xc,
		double *	yc,
		double *	zc,
		double *	rc,
		int *	ier
	)

Definition at line 9534 of file util_sparx.cpp.

References abs, circum_(), FALSE_, ierr, lstptr_(), nn(), swptst_(), t, and TRUE_.

{
    /* System generated locals */
    int i__1, i__2;

    /* Builtin functions */
    //double acos(double);

    /* Local variables */
    static double c__[3], t;
    static int i1, i2, i3, i4, n0, n1, n2, n3, n4;
    static double v1[3], v2[3], v3[3];
    static int lp, kt, nn, nt, nm2, kt1, kt2, kt11, kt12, kt21, kt22, lpl,
             lpn;
    static long int swp;
    static int ierr;
    extern /* Subroutine */ int circum_(double *, double *,
            double *, double *, int *);
    extern int lstptr_(int *, int *, int *, int *);
    extern long int swptst_(int *, int *, int *, int *,
            double *, double *, double *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   03/05/03 */

/*   Given a Delaunay triangulation of nodes on the surface */
/* of the unit sphere, this subroutine returns the set of */
/* triangle circumcenters corresponding to Voronoi vertices, */
/* along with the circumradii and a list of triangle indexes */
/* LISTC stored in one-to-one correspondence with LIST/LPTR */
/* entries. */

/*   A triangle circumcenter is the point (unit vector) lying */
/* at the same angular distance from the three vertices and */
/* contained in the same hemisphere as the vertices.  (Note */
/* that the negative of a circumcenter is also equidistant */
/* from the vertices.)  If the triangulation covers the sur- */
/* face, the Voronoi vertices are the circumcenters of the */
/* triangles in the Delaunay triangulation.  LPTR, LEND, and */
/* LNEW are not altered in this case. */

/*   On the other hand, if the nodes are contained in a sin- */
/* gle hemisphere, the triangulation is implicitly extended */
/* to the entire surface by adding pseudo-arcs (of length */
/* greater than 180 degrees) between boundary nodes forming */
/* pseudo-triangles whose 'circumcenters' are included in the */
/* list.  This extension to the triangulation actually con- */
/* sists of a triangulation of the set of boundary nodes in */
/* which the swap test is reversed (a non-empty circumcircle */
/* test).  The negative circumcenters are stored as the */
/* pseudo-triangle 'circumcenters'.  LISTC, LPTR, LEND, and */
/* LNEW contain a data structure corresponding to the ex- */
/* tended triangulation (Voronoi diagram), but LIST is not */
/* altered in this case.  Thus, if it is necessary to retain */
/* the original (unextended) triangulation data structure, */
/* copies of LPTR and LNEW must be saved before calling this */
/* routine. */


/* On input: */

/*       N = Number of nodes in the triangulation.  N .GE. 3. */
/*           Note that, if N = 3, there are only two Voronoi */
/*           vertices separated by 180 degrees, and the */
/*           Voronoi regions are not well defined. */

/*       NCOL = Number of columns reserved for LTRI.  This */
/*              must be at least NB-2, where NB is the number */
/*              of boundary nodes. */

/*       X,Y,Z = Arrays of length N containing the Cartesian */
/*               coordinates of the nodes (unit vectors). */

/*       LIST = int array containing the set of adjacency */
/*              lists.  Refer to Subroutine TRMESH. */

/*       LEND = Set of pointers to ends of adjacency lists. */
/*              Refer to Subroutine TRMESH. */

/* The above parameters are not altered by this routine. */

/*       LPTR = Array of pointers associated with LIST.  Re- */
/*              fer to Subroutine TRMESH. */

/*       LNEW = Pointer to the first empty location in LIST */
/*              and LPTR (list length plus one). */

/*       LTRI = int work space array dimensioned 6 by */
/*              NCOL, or unused dummy parameter if NB = 0. */

/*       LISTC = int array of length at least 3*NT, where */
/*               NT = 2*N-4 is the number of triangles in the */
/*               triangulation (after extending it to cover */
/*               the entire surface if necessary). */

/*       XC,YC,ZC,RC = Arrays of length NT = 2*N-4. */

/* On output: */

/*       LPTR = Array of pointers associated with LISTC: */
/*              updated for the addition of pseudo-triangles */
/*              if the original triangulation contains */
/*              boundary nodes (NB > 0). */

/*       LNEW = Pointer to the first empty location in LISTC */
/*              and LPTR (list length plus one).  LNEW is not */
/*              altered if NB = 0. */

/*       LTRI = Triangle list whose first NB-2 columns con- */
/*              tain the indexes of a clockwise-ordered */
/*              sequence of vertices (first three rows) */
/*              followed by the LTRI column indexes of the */
/*              triangles opposite the vertices (or 0 */
/*              denoting the exterior region) in the last */
/*              three rows.  This array is not generally of */
/*              any use. */

/*       LISTC = Array containing triangle indexes (indexes */
/*               to XC, YC, ZC, and RC) stored in 1-1 corres- */
/*               pondence with LIST/LPTR entries (or entries */
/*               that would be stored in LIST for the */
/*               extended triangulation):  the index of tri- */
/*               angle (N1,N2,N3) is stored in LISTC(K), */
/*               LISTC(L), and LISTC(M), where LIST(K), */
/*               LIST(L), and LIST(M) are the indexes of N2 */
/*               as a neighbor of N1, N3 as a neighbor of N2, */
/*               and N1 as a neighbor of N3.  The Voronoi */
/*               region associated with a node is defined by */
/*               the CCW-ordered sequence of circumcenters in */
/*               one-to-one correspondence with its adjacency */
/*               list (in the extended triangulation). */

/*       NB = Number of boundary nodes unless IER = 1. */

/*       XC,YC,ZC = Arrays containing the Cartesian coordi- */
/*                  nates of the triangle circumcenters */
/*                  (Voronoi vertices).  XC(I)**2 + YC(I)**2 */
/*                  + ZC(I)**2 = 1.  The first NB-2 entries */
/*                  correspond to pseudo-triangles if NB > 0. */

/*       RC = Array containing circumradii (the arc lengths */
/*            or angles between the circumcenters and associ- */
/*            ated triangle vertices) in 1-1 correspondence */
/*            with circumcenters. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if N < 3. */
/*             IER = 2 if NCOL < NB-2. */
/*             IER = 3 if a triangle is degenerate (has ver- */
/*                     tices lying on a common geodesic). */

/* Modules required by CRLIST:  CIRCUM, LSTPTR, SWPTST */

/* Intrinsic functions called by CRLIST:  ABS, ACOS */

/* *********************************************************** */


/* Local parameters: */

/* C =         Circumcenter returned by Subroutine CIRCUM */
/* I1,I2,I3 =  Permutation of (1,2,3):  LTRI row indexes */
/* I4 =        LTRI row index in the range 1 to 3 */
/* IERR =      Error flag for calls to CIRCUM */
/* KT =        Triangle index */
/* KT1,KT2 =   Indexes of a pair of adjacent pseudo-triangles */
/* KT11,KT12 = Indexes of the pseudo-triangles opposite N1 */
/*               and N2 as vertices of KT1 */
/* KT21,KT22 = Indexes of the pseudo-triangles opposite N1 */
/*               and N2 as vertices of KT2 */
/* LP,LPN =    LIST pointers */
/* LPL =       LIST pointer of the last neighbor of N1 */
/* N0 =        Index of the first boundary node (initial */
/*               value of N1) in the loop on boundary nodes */
/*               used to store the pseudo-triangle indexes */
/*               in LISTC */
/* N1,N2,N3 =  Nodal indexes defining a triangle (CCW order) */
/*               or pseudo-triangle (clockwise order) */
/* N4 =        Index of the node opposite N2 -> N1 */
/* NM2 =       N-2 */
/* NN =        Local copy of N */
/* NT =        Number of pseudo-triangles:  NB-2 */
/* SWP =       long int variable set to TRUE in each optimiza- */
/*               tion loop (loop on pseudo-arcs) iff a swap */
/*               is performed */
/* V1,V2,V3 =  Vertices of triangle KT = (N1,N2,N3) sent to */
/*               Subroutine CIRCUM */

    /* Parameter adjustments */
    --lend;
    --z__;
    --y;
    --x;
    ltri -= 7;
    --list;
    --lptr;
    --listc;
    --xc;
    --yc;
    --zc;
    --rc;

    /* Function Body */
    nn = *n;
    *nb = 0;
    nt = 0;
    if (nn < 3) {
        goto L21;
    }

/* Search for a boundary node N1. */

    i__1 = nn;
    for (n1 = 1; n1 <= i__1; ++n1) {
        lp = lend[n1];
        if (list[lp] < 0) {
            goto L2;
        }
/* L1: */
    }

/* The triangulation already covers the sphere. */

    goto L9;

/* There are NB .GE. 3 boundary nodes.  Add NB-2 pseudo- */
/*   triangles (N1,N2,N3) by connecting N3 to the NB-3 */
/*   boundary nodes to which it is not already adjacent. */

/*   Set N3 and N2 to the first and last neighbors, */
/*     respectively, of N1. */

L2:
    n2 = -list[lp];
    lp = lptr[lp];
    n3 = list[lp];

/*   Loop on boundary arcs N1 -> N2 in clockwise order, */
/*     storing triangles (N1,N2,N3) in column NT of LTRI */
/*     along with the indexes of the triangles opposite */
/*     the vertices. */

L3:
    ++nt;
    if (nt <= *ncol) {
        ltri[nt * 6 + 1] = n1;
        ltri[nt * 6 + 2] = n2;
        ltri[nt * 6 + 3] = n3;
        ltri[nt * 6 + 4] = nt + 1;
        ltri[nt * 6 + 5] = nt - 1;
        ltri[nt * 6 + 6] = 0;
    }
    n1 = n2;
    lp = lend[n1];
    n2 = -list[lp];
    if (n2 != n3) {
        goto L3;
    }

    *nb = nt + 2;
    if (*ncol < nt) {
        goto L22;
    }
    ltri[nt * 6 + 4] = 0;
    if (nt == 1) {
        goto L7;
    }

/* Optimize the exterior triangulation (set of pseudo- */
/*   triangles) by applying swaps to the pseudo-arcs N1-N2 */
/*   (pairs of adjacent pseudo-triangles KT1 and KT2 > KT1). */
/*   The loop on pseudo-arcs is repeated until no swaps are */
/*   performed. */

L4:
    swp = FALSE_;
    i__1 = nt - 1;
    for (kt1 = 1; kt1 <= i__1; ++kt1) {
        for (i3 = 1; i3 <= 3; ++i3) {
            kt2 = ltri[i3 + 3 + kt1 * 6];
            if (kt2 <= kt1) {
                goto L5;
            }

/*   The LTRI row indexes (I1,I2,I3) of triangle KT1 = */
/*     (N1,N2,N3) are a cyclical permutation of (1,2,3). */

            if (i3 == 1) {
                i1 = 2;
                i2 = 3;
            } else if (i3 == 2) {
                i1 = 3;
                i2 = 1;
            } else {
                i1 = 1;
                i2 = 2;
            }
            n1 = ltri[i1 + kt1 * 6];
            n2 = ltri[i2 + kt1 * 6];
            n3 = ltri[i3 + kt1 * 6];

/*   KT2 = (N2,N1,N4) for N4 = LTRI(I,KT2), where */
/*     LTRI(I+3,KT2) = KT1. */

            if (ltri[kt2 * 6 + 4] == kt1) {
                i4 = 1;
            } else if (ltri[kt2 * 6 + 5] == kt1) {
                i4 = 2;
            } else {
                i4 = 3;
            }
            n4 = ltri[i4 + kt2 * 6];

/*   The empty circumcircle test is reversed for the pseudo- */
/*     triangles.  The reversal is implicit in the clockwise */
/*     ordering of the vertices. */

            if (! swptst_(&n1, &n2, &n3, &n4, &x[1], &y[1], &z__[1])) {
                goto L5;
            }

/*   Swap arc N1-N2 for N3-N4.  KTij is the triangle opposite */
/*     Nj as a vertex of KTi. */

            swp = TRUE_;
            kt11 = ltri[i1 + 3 + kt1 * 6];
            kt12 = ltri[i2 + 3 + kt1 * 6];
            if (i4 == 1) {
                i2 = 2;
                i1 = 3;
            } else if (i4 == 2) {
                i2 = 3;
                i1 = 1;
            } else {
                i2 = 1;
                i1 = 2;
            }
            kt21 = ltri[i1 + 3 + kt2 * 6];
            kt22 = ltri[i2 + 3 + kt2 * 6];
            ltri[kt1 * 6 + 1] = n4;
            ltri[kt1 * 6 + 2] = n3;
            ltri[kt1 * 6 + 3] = n1;
            ltri[kt1 * 6 + 4] = kt12;
            ltri[kt1 * 6 + 5] = kt22;
            ltri[kt1 * 6 + 6] = kt2;
            ltri[kt2 * 6 + 1] = n3;
            ltri[kt2 * 6 + 2] = n4;
            ltri[kt2 * 6 + 3] = n2;
            ltri[kt2 * 6 + 4] = kt21;
            ltri[kt2 * 6 + 5] = kt11;
            ltri[kt2 * 6 + 6] = kt1;

/*   Correct the KT11 and KT22 entries that changed. */

            if (kt11 != 0) {
                i4 = 4;
                if (ltri[kt11 * 6 + 4] != kt1) {
                    i4 = 5;
                    if (ltri[kt11 * 6 + 5] != kt1) {
                        i4 = 6;
                    }
                }
                ltri[i4 + kt11 * 6] = kt2;
            }
            if (kt22 != 0) {
                i4 = 4;
                if (ltri[kt22 * 6 + 4] != kt2) {
                    i4 = 5;
                    if (ltri[kt22 * 6 + 5] != kt2) {
                        i4 = 6;
                    }
                }
                ltri[i4 + kt22 * 6] = kt1;
            }
L5:
            ;
        }
/* L6: */
    }
    if (swp) {
        goto L4;
    }

/* Compute and store the negative circumcenters and radii of */
/*   the pseudo-triangles in the first NT positions. */

L7:
    i__1 = nt;
    for (kt = 1; kt <= i__1; ++kt) {
        n1 = ltri[kt * 6 + 1];
        n2 = ltri[kt * 6 + 2];
        n3 = ltri[kt * 6 + 3];
        v1[0] = x[n1];
        v1[1] = y[n1];
        v1[2] = z__[n1];
        v2[0] = x[n2];
        v2[1] = y[n2];
        v2[2] = z__[n2];
        v3[0] = x[n3];
        v3[1] = y[n3];
        v3[2] = z__[n3];
        circum_(v2, v1, v3, c__, &ierr);
        if (ierr != 0) {
            goto L23;
        }

/*   Store the negative circumcenter and radius (computed */
/*     from <V1,C>). */

        xc[kt] = -c__[0];
        yc[kt] = -c__[1];
        zc[kt] = -c__[2];
        t = -(v1[0] * c__[0] + v1[1] * c__[1] + v1[2] * c__[2]);
        if (t < -1.) {
            t = -1.;
        }
        if (t > 1.) {
            t = 1.;
        }
        rc[kt] = acos(t);
/* L8: */
    }

/* Compute and store the circumcenters and radii of the */
/*   actual triangles in positions KT = NT+1, NT+2, ... */
/*   Also, store the triangle indexes KT in the appropriate */
/*   LISTC positions. */

L9:
    kt = nt;

/*   Loop on nodes N1. */

    nm2 = nn - 2;
    i__1 = nm2;
    for (n1 = 1; n1 <= i__1; ++n1) {
        lpl = lend[n1];
        lp = lpl;
        n3 = list[lp];

/*   Loop on adjacent neighbors N2,N3 of N1 for which N2 > N1 */
/*     and N3 > N1. */

L10:
        lp = lptr[lp];
        n2 = n3;
        n3 = (i__2 = list[lp], abs(i__2));
        if (n2 <= n1 || n3 <= n1) {
            goto L11;
        }
        ++kt;

/*   Compute the circumcenter C of triangle KT = (N1,N2,N3). */

        v1[0] = x[n1];
        v1[1] = y[n1];
        v1[2] = z__[n1];
        v2[0] = x[n2];
        v2[1] = y[n2];
        v2[2] = z__[n2];
        v3[0] = x[n3];
        v3[1] = y[n3];
        v3[2] = z__[n3];
        circum_(v1, v2, v3, c__, &ierr);
        if (ierr != 0) {
            goto L23;
        }

/*   Store the circumcenter, radius and triangle index. */

        xc[kt] = c__[0];
        yc[kt] = c__[1];
        zc[kt] = c__[2];
        t = v1[0] * c__[0] + v1[1] * c__[1] + v1[2] * c__[2];
        if (t < -1.) {
            t = -1.;
        }
        if (t > 1.) {
            t = 1.;
        }
        rc[kt] = acos(t);

/*   Store KT in LISTC(LPN), where Abs(LIST(LPN)) is the */
/*     index of N2 as a neighbor of N1, N3 as a neighbor */
/*     of N2, and N1 as a neighbor of N3. */

        lpn = lstptr_(&lpl, &n2, &list[1], &lptr[1]);
        listc[lpn] = kt;
        lpn = lstptr_(&lend[n2], &n3, &list[1], &lptr[1]);
        listc[lpn] = kt;
        lpn = lstptr_(&lend[n3], &n1, &list[1], &lptr[1]);
        listc[lpn] = kt;
L11:
        if (lp != lpl) {
            goto L10;
        }
/* L12: */
    }
    if (nt == 0) {
        goto L20;
    }

/* Store the first NT triangle indexes in LISTC. */

/*   Find a boundary triangle KT1 = (N1,N2,N3) with a */
/*     boundary arc opposite N3. */

    kt1 = 0;
L13:
    ++kt1;
    if (ltri[kt1 * 6 + 4] == 0) {
        i1 = 2;
        i2 = 3;
        i3 = 1;
        goto L14;
    } else if (ltri[kt1 * 6 + 5] == 0) {
        i1 = 3;
        i2 = 1;
        i3 = 2;
        goto L14;
    } else if (ltri[kt1 * 6 + 6] == 0) {
        i1 = 1;
        i2 = 2;
        i3 = 3;
        goto L14;
    }
    goto L13;
L14:
    n1 = ltri[i1 + kt1 * 6];
    n0 = n1;

/*   Loop on boundary nodes N1 in CCW order, storing the */
/*     indexes of the clockwise-ordered sequence of triangles */
/*     that contain N1.  The first triangle overwrites the */
/*     last neighbor position, and the remaining triangles, */
/*     if any, are appended to N1's adjacency list. */

/*   A pointer to the first neighbor of N1 is saved in LPN. */

L15:
    lp = lend[n1];
    lpn = lptr[lp];
    listc[lp] = kt1;

/*   Loop on triangles KT2 containing N1. */

L16:
    kt2 = ltri[i2 + 3 + kt1 * 6];
    if (kt2 != 0) {

/*   Append KT2 to N1's triangle list. */

        lptr[lp] = *lnew;
        lp = *lnew;
        listc[lp] = kt2;
        ++(*lnew);

/*   Set KT1 to KT2 and update (I1,I2,I3) such that */
/*     LTRI(I1,KT1) = N1. */

        kt1 = kt2;
        if (ltri[kt1 * 6 + 1] == n1) {
            i1 = 1;
            i2 = 2;
            i3 = 3;
        } else if (ltri[kt1 * 6 + 2] == n1) {
            i1 = 2;
            i2 = 3;
            i3 = 1;
        } else {
            i1 = 3;
            i2 = 1;
            i3 = 2;
        }
        goto L16;
    }

/*   Store the saved first-triangle pointer in LPTR(LP), set */
/*     N1 to the next boundary node, test for termination, */
/*     and permute the indexes:  the last triangle containing */
/*     a boundary node is the first triangle containing the */
/*     next boundary node. */

    lptr[lp] = lpn;
    n1 = ltri[i3 + kt1 * 6];
    if (n1 != n0) {
        i4 = i3;
        i3 = i2;
        i2 = i1;
        i1 = i4;
        goto L15;
    }

/* No errors encountered. */

L20:
    *ier = 0;
    return 0;

/* N < 3. */

L21:
    *ier = 1;
    return 0;

/* Insufficient space reserved for LTRI. */

L22:
    *ier = 2;
    return 0;

/* Error flag returned by CIRCUM: KT indexes a null triangle. */

L23:
    *ier = 3;
    return 0;
} /* crlist_ */

int delarc_	(	int *	n,
		int *	io1,
		int *	io2,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	lnew,
		int *	ier
	)

Definition at line 10164 of file util_sparx.cpp.

References abs, delnb_(), and lstptr_().

{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int n1, n2, n3, lp, lph, lpl;
    extern /* Subroutine */ int delnb_(int *, int *, int *,
            int *, int *, int *, int *, int *);
    extern int lstptr_(int *, int *, int *, int *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/17/96 */

/*   This subroutine deletes a boundary arc from a triangula- */
/* tion.  It may be used to remove a null triangle from the */
/* convex hull boundary.  Note, however, that if the union of */
/* triangles is rendered nonconvex, Subroutines DELNOD, EDGE, */
/* and TRFIND (and hence ADDNOD) may fail.  Also, Function */
/* NEARND should not be called following an arc deletion. */

/*   This routine is identical to the similarly named routine */
/* in TRIPACK. */


/* On input: */

/*       N = Number of nodes in the triangulation.  N .GE. 4. */

/*       IO1,IO2 = Indexes (in the range 1 to N) of a pair of */
/*                 adjacent boundary nodes defining the arc */
/*                 to be removed. */

/* The above parameters are not altered by this routine. */

/*       LIST,LPTR,LEND,LNEW = Triangulation data structure */
/*                             created by Subroutine TRMESH. */

/* On output: */

/*       LIST,LPTR,LEND,LNEW = Data structure updated with */
/*                             the removal of arc IO1-IO2 */
/*                             unless IER > 0. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if N, IO1, or IO2 is outside its valid */
/*                     range, or IO1 = IO2. */
/*             IER = 2 if IO1-IO2 is not a boundary arc. */
/*             IER = 3 if the node opposite IO1-IO2 is al- */
/*                     ready a boundary node, and thus IO1 */
/*                     or IO2 has only two neighbors or a */
/*                     deletion would result in two triangu- */
/*                     lations sharing a single node. */
/*             IER = 4 if one of the nodes is a neighbor of */
/*                     the other, but not vice versa, imply- */
/*                     ing an invalid triangulation data */
/*                     structure. */

/* Module required by DELARC:  DELNB, LSTPTR */

/* Intrinsic function called by DELARC:  ABS */

/* *********************************************************** */


/* Local parameters: */

/* LP =       LIST pointer */
/* LPH =      LIST pointer or flag returned by DELNB */
/* LPL =      Pointer to the last neighbor of N1, N2, or N3 */
/* N1,N2,N3 = Nodal indexes of a triangle such that N1->N2 */
/*              is the directed boundary edge associated */
/*              with IO1-IO2 */

    /* Parameter adjustments */
    --lend;
    --list;
    --lptr;

    /* Function Body */
    n1 = *io1;
    n2 = *io2;

/* Test for errors, and set N1->N2 to the directed boundary */
/*   edge associated with IO1-IO2:  (N1,N2,N3) is a triangle */
/*   for some N3. */

    if (*n < 4 || n1 < 1 || n1 > *n || n2 < 1 || n2 > *n || n1 == n2) {
        *ier = 1;
        return 0;
    }

    lpl = lend[n2];
    if (-list[lpl] != n1) {
        n1 = n2;
        n2 = *io1;
        lpl = lend[n2];
        if (-list[lpl] != n1) {
            *ier = 2;
            return 0;
        }
    }

/* Set N3 to the node opposite N1->N2 (the second neighbor */
/*   of N1), and test for error 3 (N3 already a boundary */
/*   node). */

    lpl = lend[n1];
    lp = lptr[lpl];
    lp = lptr[lp];
    n3 = (i__1 = list[lp], abs(i__1));
    lpl = lend[n3];
    if (list[lpl] <= 0) {
        *ier = 3;
        return 0;
    }

/* Delete N2 as a neighbor of N1, making N3 the first */
/*   neighbor, and test for error 4 (N2 not a neighbor */
/*   of N1).  Note that previously computed pointers may */
/*   no longer be valid following the call to DELNB. */

    delnb_(&n1, &n2, n, &list[1], &lptr[1], &lend[1], lnew, &lph);
    if (lph < 0) {
        *ier = 4;
        return 0;
    }

/* Delete N1 as a neighbor of N2, making N3 the new last */
/*   neighbor. */

    delnb_(&n2, &n1, n, &list[1], &lptr[1], &lend[1], lnew, &lph);

/* Make N3 a boundary node with first neighbor N2 and last */
/*   neighbor N1. */

    lp = lstptr_(&lend[n3], &n1, &list[1], &lptr[1]);
    lend[n3] = lp;
    list[lp] = -n1;

/* No errors encountered. */

    *ier = 0;
    return 0;
} /* delarc_ */

int delnb_	(	int *	n0,
		int *	nb,
		int *	n,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	lnew,
		int *	lph
	)

Definition at line 10319 of file util_sparx.cpp.

References abs, and nn().

Referenced by delarc_(), and delnod_().

{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int i__, lp, nn, lpb, lpl, lpp, lnw;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/29/98 */

/*   This subroutine deletes a neighbor NB from the adjacency */
/* list of node N0 (but N0 is not deleted from the adjacency */
/* list of NB) and, if NB is a boundary node, makes N0 a */
/* boundary node.  For pointer (LIST index) LPH to NB as a */
/* neighbor of N0, the empty LIST,LPTR location LPH is filled */
/* in with the values at LNEW-1, pointer LNEW-1 (in LPTR and */
/* possibly in LEND) is changed to LPH, and LNEW is decremen- */
/* ted.  This requires a search of LEND and LPTR entailing an */
/* expected operation count of O(N). */

/*   This routine is identical to the similarly named routine */
/* in TRIPACK. */


/* On input: */

/*       N0,NB = Indexes, in the range 1 to N, of a pair of */
/*               nodes such that NB is a neighbor of N0. */
/*               (N0 need not be a neighbor of NB.) */

/*       N = Number of nodes in the triangulation.  N .GE. 3. */

/* The above parameters are not altered by this routine. */

/*       LIST,LPTR,LEND,LNEW = Data structure defining the */
/*                             triangulation. */

/* On output: */

/*       LIST,LPTR,LEND,LNEW = Data structure updated with */
/*                             the removal of NB from the ad- */
/*                             jacency list of N0 unless */
/*                             LPH < 0. */

/*       LPH = List pointer to the hole (NB as a neighbor of */
/*             N0) filled in by the values at LNEW-1 or error */
/*             indicator: */
/*             LPH > 0 if no errors were encountered. */
/*             LPH = -1 if N0, NB, or N is outside its valid */
/*                      range. */
/*             LPH = -2 if NB is not a neighbor of N0. */

/* Modules required by DELNB:  None */

/* Intrinsic function called by DELNB:  ABS */

/* *********************************************************** */


/* Local parameters: */

/* I =   DO-loop index */
/* LNW = LNEW-1 (output value of LNEW) */
/* LP =  LIST pointer of the last neighbor of NB */
/* LPB = Pointer to NB as a neighbor of N0 */
/* LPL = Pointer to the last neighbor of N0 */
/* LPP = Pointer to the neighbor of N0 that precedes NB */
/* NN =  Local copy of N */

    /* Parameter adjustments */
    --lend;
    --list;
    --lptr;

    /* Function Body */
    nn = *n;

/* Test for error 1. */

    if (*n0 < 1 || *n0 > nn || *nb < 1 || *nb > nn || nn < 3) {
        *lph = -1;
        return 0;
    }

/*   Find pointers to neighbors of N0: */

/*     LPL points to the last neighbor, */
/*     LPP points to the neighbor NP preceding NB, and */
/*     LPB points to NB. */

    lpl = lend[*n0];
    lpp = lpl;
    lpb = lptr[lpp];
L1:
    if (list[lpb] == *nb) {
        goto L2;
    }
    lpp = lpb;
    lpb = lptr[lpp];
    if (lpb != lpl) {
        goto L1;
    }

/*   Test for error 2 (NB not found). */

    if ((i__1 = list[lpb], abs(i__1)) != *nb) {
        *lph = -2;
        return 0;
    }

/*   NB is the last neighbor of N0.  Make NP the new last */
/*     neighbor and, if NB is a boundary node, then make N0 */
/*     a boundary node. */

    lend[*n0] = lpp;
    lp = lend[*nb];
    if (list[lp] < 0) {
        list[lpp] = -list[lpp];
    }
    goto L3;

/*   NB is not the last neighbor of N0.  If NB is a boundary */
/*     node and N0 is not, then make N0 a boundary node with */
/*     last neighbor NP. */

L2:
    lp = lend[*nb];
    if (list[lp] < 0 && list[lpl] > 0) {
        lend[*n0] = lpp;
        list[lpp] = -list[lpp];
    }

/*   Update LPTR so that the neighbor following NB now fol- */
/*     lows NP, and fill in the hole at location LPB. */

L3:
    lptr[lpp] = lptr[lpb];
    lnw = *lnew - 1;
    list[lpb] = list[lnw];
    lptr[lpb] = lptr[lnw];
    for (i__ = nn; i__ >= 1; --i__) {
        if (lend[i__] == lnw) {
            lend[i__] = lpb;
            goto L5;
        }
/* L4: */
    }

L5:
    i__1 = lnw - 1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        if (lptr[i__] == lnw) {
            lptr[i__] = lpb;
        }
/* L6: */
    }

/* No errors encountered. */

    *lnew = lnw;
    *lph = lpb;
    return 0;
} /* delnb_ */

int delnod_	(	int *	k,
		int *	n,
		double *	x,
		double *	y,
		double *	z__,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	lnew,
		int *	lwk,
		int *	iwk,
		int *	ier
	)

Definition at line 10492 of file util_sparx.cpp.

References abs, delnb_(), FALSE_, ierr, left_(), lstptr_(), nbcnt_(), nn(), optim_(), swap_(), and TRUE_.

{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int i__, j, n1, n2;
    static double x1, x2, y1, y2, z1, z2;
    static int nl, lp, nn, nr;
    static double xl, yl, zl, xr, yr, zr;
    static int nnb, lp21, lpf, lph, lpl, lpn, iwl, nit, lnw, lpl2;
    extern long int left_(double *, double *, double *, double
            *, double *, double *, double *, double *,
            double *);
    static long int bdry;
    static int ierr, lwkl;
    extern /* Subroutine */ int swap_(int *, int *, int *,
            int *, int *, int *, int *, int *), delnb_(
            int *, int *, int *, int *, int *, int *,
            int *, int *);
    extern int nbcnt_(int *, int *);
    extern /* Subroutine */ int optim_(double *, double *, double
            *, int *, int *, int *, int *, int *, int
            *, int *);
    static int nfrst;
    extern int lstptr_(int *, int *, int *, int *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   11/30/99 */

/*   This subroutine deletes node K (along with all arcs */
/* incident on node K) from a triangulation of N nodes on the */
/* unit sphere, and inserts arcs as necessary to produce a */
/* triangulation of the remaining N-1 nodes.  If a Delaunay */
/* triangulation is input, a Delaunay triangulation will */
/* result, and thus, DELNOD reverses the effect of a call to */
/* Subroutine ADDNOD. */


/* On input: */

/*       K = Index (for X, Y, and Z) of the node to be */
/*           deleted.  1 .LE. K .LE. N. */

/* K is not altered by this routine. */

/*       N = Number of nodes in the triangulation on input. */
/*           N .GE. 4.  Note that N will be decremented */
/*           following the deletion. */

/*       X,Y,Z = Arrays of length N containing the Cartesian */
/*               coordinates of the nodes in the triangula- */
/*               tion. */

/*       LIST,LPTR,LEND,LNEW = Data structure defining the */
/*                             triangulation.  Refer to Sub- */
/*                             routine TRMESH. */

/*       LWK = Number of columns reserved for IWK.  LWK must */
/*             be at least NNB-3, where NNB is the number of */
/*             neighbors of node K, including an extra */
/*             pseudo-node if K is a boundary node. */

/*       IWK = int work array dimensioned 2 by LWK (or */
/*             array of length .GE. 2*LWK). */

/* On output: */

/*       N = Number of nodes in the triangulation on output. */
/*           The input value is decremented unless 1 .LE. IER */
/*           .LE. 4. */

/*       X,Y,Z = Updated arrays containing nodal coordinates */
/*               (with elements K+1,...,N+1 shifted up one */
/*               position, thus overwriting element K) unless */
/*               1 .LE. IER .LE. 4. */

/*       LIST,LPTR,LEND,LNEW = Updated triangulation data */
/*                             structure reflecting the dele- */
/*                             tion unless 1 .LE. IER .LE. 4. */
/*                             Note that the data structure */
/*                             may have been altered if IER > */
/*                             3. */

/*       LWK = Number of IWK columns required unless IER = 1 */
/*             or IER = 3. */

/*       IWK = Indexes of the endpoints of the new arcs added */
/*             unless LWK = 0 or 1 .LE. IER .LE. 4.  (Arcs */
/*             are associated with columns, or pairs of */
/*             adjacent elements if IWK is declared as a */
/*             singly-subscripted array.) */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if K or N is outside its valid range */
/*                     or LWK < 0 on input. */
/*             IER = 2 if more space is required in IWK. */
/*                     Refer to LWK. */
/*             IER = 3 if the triangulation data structure is */
/*                     invalid on input. */
/*             IER = 4 if K indexes an interior node with */
/*                     four or more neighbors, none of which */
/*                     can be swapped out due to collineari- */
/*                     ty, and K cannot therefore be deleted. */
/*             IER = 5 if an error flag (other than IER = 1) */
/*                     was returned by OPTIM.  An error */
/*                     message is written to the standard */
/*                     output unit in this case. */
/*             IER = 6 if error flag 1 was returned by OPTIM. */
/*                     This is not necessarily an error, but */
/*                     the arcs may not be optimal. */

/*   Note that the deletion may result in all remaining nodes */
/* being collinear.  This situation is not flagged. */

/* Modules required by DELNOD:  DELNB, LEFT, LSTPTR, NBCNT, */
/*                                OPTIM, SWAP, SWPTST */

/* Intrinsic function called by DELNOD:  ABS */

/* *********************************************************** */


/* Local parameters: */

/* BDRY =    long int variable with value TRUE iff N1 is a */
/*             boundary node */
/* I,J =     DO-loop indexes */
/* IERR =    Error flag returned by OPTIM */
/* IWL =     Number of IWK columns containing arcs */
/* LNW =     Local copy of LNEW */
/* LP =      LIST pointer */
/* LP21 =    LIST pointer returned by SWAP */
/* LPF,LPL = Pointers to the first and last neighbors of N1 */
/* LPH =     Pointer (or flag) returned by DELNB */
/* LPL2 =    Pointer to the last neighbor of N2 */
/* LPN =     Pointer to a neighbor of N1 */
/* LWKL =    Input value of LWK */
/* N1 =      Local copy of K */
/* N2 =      Neighbor of N1 */
/* NFRST =   First neighbor of N1:  LIST(LPF) */
/* NIT =     Number of iterations in OPTIM */
/* NR,NL =   Neighbors of N1 preceding (to the right of) and */
/*             following (to the left of) N2, respectively */
/* NN =      Number of nodes in the triangulation */
/* NNB =     Number of neighbors of N1 (including a pseudo- */
/*             node representing the boundary if N1 is a */
/*             boundary node) */
/* X1,Y1,Z1 = Coordinates of N1 */
/* X2,Y2,Z2 = Coordinates of N2 */
/* XL,YL,ZL = Coordinates of NL */
/* XR,YR,ZR = Coordinates of NR */


/* Set N1 to K and NNB to the number of neighbors of N1 (plus */
/*   one if N1 is a boundary node), and test for errors.  LPF */
/*   and LPL are LIST indexes of the first and last neighbors */
/*   of N1, IWL is the number of IWK columns containing arcs, */
/*   and BDRY is TRUE iff N1 is a boundary node. */

    /* Parameter adjustments */
    iwk -= 3;
    --lend;
    --lptr;
    --list;
    --z__;
    --y;
    --x;

    /* Function Body */
    n1 = *k;
    nn = *n;
    if (n1 < 1 || n1 > nn || nn < 4 || *lwk < 0) {
        goto L21;
    }
    lpl = lend[n1];
    lpf = lptr[lpl];
    nnb = nbcnt_(&lpl, &lptr[1]);
    bdry = list[lpl] < 0;
    if (bdry) {
        ++nnb;
    }
    if (nnb < 3) {
        goto L23;
    }
    lwkl = *lwk;
    *lwk = nnb - 3;
    if (lwkl < *lwk) {
        goto L22;
    }
    iwl = 0;
    if (nnb == 3) {
        goto L3;
    }

/* Initialize for loop on arcs N1-N2 for neighbors N2 of N1, */
/*   beginning with the second neighbor.  NR and NL are the */
/*   neighbors preceding and following N2, respectively, and */
/*   LP indexes NL.  The loop is exited when all possible */
/*   swaps have been applied to arcs incident on N1. */

    x1 = x[n1];
    y1 = y[n1];
    z1 = z__[n1];
    nfrst = list[lpf];
    nr = nfrst;
    xr = x[nr];
    yr = y[nr];
    zr = z__[nr];
    lp = lptr[lpf];
    n2 = list[lp];
    x2 = x[n2];
    y2 = y[n2];
    z2 = z__[n2];
    lp = lptr[lp];

/* Top of loop:  set NL to the neighbor following N2. */

L1:
    nl = (i__1 = list[lp], abs(i__1));
    if (nl == nfrst && bdry) {
        goto L3;
    }
    xl = x[nl];
    yl = y[nl];
    zl = z__[nl];

/*   Test for a convex quadrilateral.  To avoid an incorrect */
/*     test caused by collinearity, use the fact that if N1 */
/*     is a boundary node, then N1 LEFT NR->NL and if N2 is */
/*     a boundary node, then N2 LEFT NL->NR. */

    lpl2 = lend[n2];
    if (! ((bdry || left_(&xr, &yr, &zr, &xl, &yl, &zl, &x1, &y1, &z1)) && (
            list[lpl2] < 0 || left_(&xl, &yl, &zl, &xr, &yr, &zr, &x2, &y2, &
            z2)))) {

/*   Nonconvex quadrilateral -- no swap is possible. */

        nr = n2;
        xr = x2;
        yr = y2;
        zr = z2;
        goto L2;
    }

/*   The quadrilateral defined by adjacent triangles */
/*     (N1,N2,NL) and (N2,N1,NR) is convex.  Swap in */
/*     NL-NR and store it in IWK unless NL and NR are */
/*     already adjacent, in which case the swap is not */
/*     possible.  Indexes larger than N1 must be decremented */
/*     since N1 will be deleted from X, Y, and Z. */

    swap_(&nl, &nr, &n1, &n2, &list[1], &lptr[1], &lend[1], &lp21);
    if (lp21 == 0) {
        nr = n2;
        xr = x2;
        yr = y2;
        zr = z2;
        goto L2;
    }
    ++iwl;
    if (nl <= n1) {
        iwk[(iwl << 1) + 1] = nl;
    } else {
        iwk[(iwl << 1) + 1] = nl - 1;
    }
    if (nr <= n1) {
        iwk[(iwl << 1) + 2] = nr;
    } else {
        iwk[(iwl << 1) + 2] = nr - 1;
    }

/*   Recompute the LIST indexes and NFRST, and decrement NNB. */

    lpl = lend[n1];
    --nnb;
    if (nnb == 3) {
        goto L3;
    }
    lpf = lptr[lpl];
    nfrst = list[lpf];
    lp = lstptr_(&lpl, &nl, &list[1], &lptr[1]);
    if (nr == nfrst) {
        goto L2;
    }

/*   NR is not the first neighbor of N1. */
/*     Back up and test N1-NR for a swap again:  Set N2 to */
/*     NR and NR to the previous neighbor of N1 -- the */
/*     neighbor of NR which follows N1.  LP21 points to NL */
/*     as a neighbor of NR. */

    n2 = nr;
    x2 = xr;
    y2 = yr;
    z2 = zr;
    lp21 = lptr[lp21];
    lp21 = lptr[lp21];
    nr = (i__1 = list[lp21], abs(i__1));
    xr = x[nr];
    yr = y[nr];
    zr = z__[nr];
    goto L1;

/*   Bottom of loop -- test for termination of loop. */

L2:
    if (n2 == nfrst) {
        goto L3;
    }
    n2 = nl;
    x2 = xl;
    y2 = yl;
    z2 = zl;
    lp = lptr[lp];
    goto L1;

/* Delete N1 and all its incident arcs.  If N1 is an interior */
/*   node and either NNB > 3 or NNB = 3 and N2 LEFT NR->NL, */
/*   then N1 must be separated from its neighbors by a plane */
/*   containing the origin -- its removal reverses the effect */
/*   of a call to COVSPH, and all its neighbors become */
/*   boundary nodes.  This is achieved by treating it as if */
/*   it were a boundary node (setting BDRY to TRUE, changing */
/*   a sign in LIST, and incrementing NNB). */

L3:
    if (! bdry) {
        if (nnb > 3) {
            bdry = TRUE_;
        } else {
            lpf = lptr[lpl];
            nr = list[lpf];
            lp = lptr[lpf];
            n2 = list[lp];
            nl = list[lpl];
            bdry = left_(&x[nr], &y[nr], &z__[nr], &x[nl], &y[nl], &z__[nl], &
                    x[n2], &y[n2], &z__[n2]);
        }
        if (bdry) {

/*   IF a boundary node already exists, then N1 and its */
/*     neighbors cannot be converted to boundary nodes. */
/*     (They must be collinear.)  This is a problem if */
/*     NNB > 3. */

            i__1 = nn;
            for (i__ = 1; i__ <= i__1; ++i__) {
                if (list[lend[i__]] < 0) {
                    bdry = FALSE_;
                    goto L5;
                }
/* L4: */
            }
            list[lpl] = -list[lpl];
            ++nnb;
        }
    }
L5:
    if (! bdry && nnb > 3) {
        goto L24;
    }

/* Initialize for loop on neighbors.  LPL points to the last */
/*   neighbor of N1.  LNEW is stored in local variable LNW. */

    lp = lpl;
    lnw = *lnew;

/* Loop on neighbors N2 of N1, beginning with the first. */

L6:
    lp = lptr[lp];
    n2 = (i__1 = list[lp], abs(i__1));
    delnb_(&n2, &n1, n, &list[1], &lptr[1], &lend[1], &lnw, &lph);
    if (lph < 0) {
        goto L23;
    }

/*   LP and LPL may require alteration. */

    if (lpl == lnw) {
        lpl = lph;
    }
    if (lp == lnw) {
        lp = lph;
    }
    if (lp != lpl) {
        goto L6;
    }

/* Delete N1 from X, Y, Z, and LEND, and remove its adjacency */
/*   list from LIST and LPTR.  LIST entries (nodal indexes) */
/*   which are larger than N1 must be decremented. */

    --nn;
    if (n1 > nn) {
        goto L9;
    }
    i__1 = nn;
    for (i__ = n1; i__ <= i__1; ++i__) {
        x[i__] = x[i__ + 1];
        y[i__] = y[i__ + 1];
        z__[i__] = z__[i__ + 1];
        lend[i__] = lend[i__ + 1];
/* L7: */
    }

    i__1 = lnw - 1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        if (list[i__] > n1) {
            --list[i__];
        }
        if (list[i__] < -n1) {
            ++list[i__];
        }
/* L8: */
    }

/*   For LPN = first to last neighbors of N1, delete the */
/*     preceding neighbor (indexed by LP). */

/*   Each empty LIST,LPTR location LP is filled in with the */
/*     values at LNW-1, and LNW is decremented.  All pointers */
/*     (including those in LPTR and LEND) with value LNW-1 */
/*     must be changed to LP. */

/*  LPL points to the last neighbor of N1. */

L9:
    if (bdry) {
        --nnb;
    }
    lpn = lpl;
    i__1 = nnb;
    for (j = 1; j <= i__1; ++j) {
        --lnw;
        lp = lpn;
        lpn = lptr[lp];
        list[lp] = list[lnw];
        lptr[lp] = lptr[lnw];
        if (lptr[lpn] == lnw) {
            lptr[lpn] = lp;
        }
        if (lpn == lnw) {
            lpn = lp;
        }
        for (i__ = nn; i__ >= 1; --i__) {
            if (lend[i__] == lnw) {
                lend[i__] = lp;
                goto L11;
            }
/* L10: */
        }

L11:
        for (i__ = lnw - 1; i__ >= 1; --i__) {
            if (lptr[i__] == lnw) {
                lptr[i__] = lp;
            }
/* L12: */
        }
/* L13: */
    }

/* Update N and LNEW, and optimize the patch of triangles */
/*   containing K (on input) by applying swaps to the arcs */
/*   in IWK. */

    *n = nn;
    *lnew = lnw;
    if (iwl > 0) {
        nit = iwl << 2;
        optim_(&x[1], &y[1], &z__[1], &iwl, &list[1], &lptr[1], &lend[1], &
                nit, &iwk[3], &ierr);
        if (ierr != 0 && ierr != 1) {
            goto L25;
        }
        if (ierr == 1) {
            goto L26;
        }
    }

/* Successful termination. */

    *ier = 0;
    return 0;

/* Invalid input parameter. */

L21:
    *ier = 1;
    return 0;

/* Insufficient space reserved for IWK. */

L22:
    *ier = 2;
    return 0;

/* Invalid triangulation data structure.  NNB < 3 on input or */
/*   N2 is a neighbor of N1 but N1 is not a neighbor of N2. */

L23:
    *ier = 3;
    return 0;

/* N1 is interior but NNB could not be reduced to 3. */

L24:
    *ier = 4;
    return 0;

/* Error flag (other than 1) returned by OPTIM. */

L25:
    *ier = 5;
/*      WRITE (*,100) NIT, IERR */
/*  100 FORMAT (//5X,'*** Error in OPTIM (called from ', */
/*     .        'DELNOD):  NIT = ',I4,', IER = ',I1,' ***'/) */
    return 0;

/* Error flag 1 returned by OPTIM. */

L26:
    *ier = 6;
    return 0;
} /* delnod_ */

int drwarc_	(	int *	,
		double *	p,
		double *	q,
		double *	tol,
		int *	nseg
	)

Definition at line 11032 of file util_sparx.cpp.

References abs, and sqrt().

Referenced by trplot_(), and vrplot_().

{
    /* System generated locals */
    int i__1;
    double d__1;

    /* Builtin functions */
    //double sqrt(double);

    /* Local variables */
    static int i__, k;
    static double s, p1[3], p2[3], u1, u2, v1, v2;
    static int na;
    static double dp[3], du, dv, pm[3], um, vm, err, enrm;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   03/04/03 */

/*   Given unit vectors P and Q corresponding to northern */
/* hemisphere points (with positive third components), this */
/* subroutine draws a polygonal line which approximates the */
/* projection of arc P-Q onto the plane containing the */
/* equator. */

/*   The line segment is drawn by writing a sequence of */
/* 'moveto' and 'lineto' Postscript commands to unit LUN.  It */
/* is assumed that an open file is attached to the unit, */
/* header comments have been written to the file, a window- */
/* to-viewport mapping has been established, etc. */

/* On input: */

/*       LUN = long int unit number in the range 0 to 99. */

/*       P,Q = Arrays of length 3 containing the endpoints of */
/*             the arc to be drawn. */

/*       TOL = Maximum distance in world coordinates between */
/*             the projected arc and polygonal line. */

/* Input parameters are not altered by this routine. */

/* On output: */

/*       NSEG = Number of line segments in the polygonal */
/*              approximation to the projected arc.  This is */
/*              a decreasing function of TOL.  NSEG = 0 and */
/*              no drawing is performed if P = Q or P = -Q */
/*              or an error is encountered in writing to unit */
/*              LUN. */

/* STRIPACK modules required by DRWARC:  None */

/* Intrinsic functions called by DRWARC:  ABS, DBLE, SQRT */

/* *********************************************************** */


/* Local parameters: */

/* DP =    (Q-P)/NSEG */
/* DU,DV = Components of the projection Q'-P' of arc P->Q */
/*           onto the projection plane */
/* ENRM =  Euclidean norm (or squared norm) of Q'-P' or PM */
/* ERR =   Orthogonal distance from the projected midpoint */
/*           PM' to the line defined by P' and Q': */
/*           |Q'-P' X PM'-P'|/|Q'-P'| */
/* I,K =   DO-loop indexes */
/* NA =    Number of arcs (segments) in the partition of P-Q */
/* P1,P2 = Pairs of adjacent points in a uniform partition of */
/*           arc P-Q into NSEG segments; obtained by normal- */
/*           izing PM values */
/* PM =    Midpoint of arc P-Q or a point P + k*DP in a */
/*           uniform partition of the line segment P-Q into */
/*           NSEG segments */
/* S =     Scale factor 1/NA */
/* U1,V1 = Components of P' */
/* U2,V2 = Components of Q' */
/* UM,VM = Components of the midpoint PM' */


/* Compute the midpoint PM of arc P-Q. */

    /* Parameter adjustments */
    --q;
    --p;

    /* Function Body */
    enrm = 0.;
    for (i__ = 1; i__ <= 3; ++i__) {
        pm[i__ - 1] = p[i__] + q[i__];
        enrm += pm[i__ - 1] * pm[i__ - 1];
/* L1: */
    }
    if (enrm == 0.) {
        goto L5;
    }
    enrm = sqrt(enrm);
    pm[0] /= enrm;
    pm[1] /= enrm;
    pm[2] /= enrm;

/* Project P, Q, and PM to P' = (U1,V1), Q' = (U2,V2), and */
/*   PM' = (UM,VM), respectively. */

    u1 = p[1];
    v1 = p[2];
    u2 = q[1];
    v2 = q[2];
    um = pm[0];
    vm = pm[1];

/* Compute the orthogonal distance ERR from PM' to the line */
/*   defined by P' and Q'.  This is the maximum deviation */
/*   between the projected arc and the line segment.  It is */
/*   undefined if P' = Q'. */

    du = u2 - u1;
    dv = v2 - v1;
    enrm = du * du + dv * dv;
    if (enrm == 0.) {
        goto L5;
    }
    err = (d__1 = du * (vm - v1) - (um - u1) * dv, abs(d__1)) / sqrt(enrm);

/* Compute the number of arcs into which P-Q will be parti- */
/*   tioned (the number of line segments to be drawn): */
/*   NA = ERR/TOL. */

    na = (int) (err / *tol + 1.);

/* Initialize for loop on arcs P1-P2, where the intermediate */
/*   points are obtained by normalizing PM = P + k*DP for */
/*   DP = (Q-P)/NA */

    s = 1. / (double) na;
    for (i__ = 1; i__ <= 3; ++i__) {
        dp[i__ - 1] = s * (q[i__] - p[i__]);
        pm[i__ - 1] = p[i__];
        p1[i__ - 1] = p[i__];
/* L2: */
    }

/* Loop on arcs P1-P2, drawing the line segments associated */
/*   with the projected endpoints. */

    i__1 = na - 1;
    for (k = 1; k <= i__1; ++k) {
        enrm = 0.;
        for (i__ = 1; i__ <= 3; ++i__) {
            pm[i__ - 1] += dp[i__ - 1];
            enrm += pm[i__ - 1] * pm[i__ - 1];
/* L3: */
        }
        if (enrm == 0.) {
            goto L5;
        }
        enrm = sqrt(enrm);
        p2[0] = pm[0] / enrm;
        p2[1] = pm[1] / enrm;
        p2[2] = pm[2] / enrm;
/*        WRITE (LUN,100,ERR=5) P1(1), P1(2), P2(1), P2(2) */
/*  100   FORMAT (2F12.6,' moveto',2F12.6,' lineto') */
        p1[0] = p2[0];
        p1[1] = p2[1];
        p1[2] = p2[2];
/* L4: */
    }
/*      WRITE (LUN,100,ERR=5) P1(1), P1(2), Q(1), Q(2) */

/* No error encountered. */

    *nseg = na;
    return 0;

/* Invalid input value of P or Q. */

L5:
    *nseg = 0;
    return 0;
} /* drwarc_ */

int edge_	(	int *	in1,
		int *	in2,
		double *	x,
		double *	y,
		double *	z__,
		int *	lwk,
		int *	iwk,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	ier
	)

Definition at line 11222 of file util_sparx.cpp.

References abs, ierr, left_(), optim_(), and swap_().

{
    /* System generated locals */
    int i__1;

    /* Local variables */
    static int i__, n0, n1, n2;
    static double x0, x1, x2, y0, y1, y2, z0, z1, z2;
    static int nl, lp, nr;
    static double dp12;
    static int lp21, iwc, iwf, lft, lpl, iwl, nit;
    static double dp1l, dp2l, dp1r, dp2r;
    extern long int left_(double *, double *, double *, double
            *, double *, double *, double *, double *,
            double *);
    static int ierr;
    extern /* Subroutine */ int swap_(int *, int *, int *,
            int *, int *, int *, int *, int *);
    static int next, iwcp1, n1lst, iwend;
    extern /* Subroutine */ int optim_(double *, double *, double
            *, int *, int *, int *, int *, int *, int
            *, int *);
    static int n1frst;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/30/98 */

/*   Given a triangulation of N nodes and a pair of nodal */
/* indexes IN1 and IN2, this routine swaps arcs as necessary */
/* to force IN1 and IN2 to be adjacent.  Only arcs which */
/* intersect IN1-IN2 are swapped out.  If a Delaunay triangu- */
/* lation is input, the resulting triangulation is as close */
/* as possible to a Delaunay triangulation in the sense that */
/* all arcs other than IN1-IN2 are locally optimal. */

/*   A sequence of calls to EDGE may be used to force the */
/* presence of a set of edges defining the boundary of a non- */
/* convex and/or multiply connected region, or to introduce */
/* barriers into the triangulation.  Note that Subroutine */
/* GETNP will not necessarily return closest nodes if the */
/* triangulation has been constrained by a call to EDGE. */
/* However, this is appropriate in some applications, such */
/* as triangle-based interpolation on a nonconvex domain. */


/* On input: */

/*       IN1,IN2 = Indexes (of X, Y, and Z) in the range 1 to */
/*                 N defining a pair of nodes to be connected */
/*                 by an arc. */

/*       X,Y,Z = Arrays of length N containing the Cartesian */
/*               coordinates of the nodes. */

/* The above parameters are not altered by this routine. */

/*       LWK = Number of columns reserved for IWK.  This must */
/*             be at least NI -- the number of arcs that */
/*             intersect IN1-IN2.  (NI is bounded by N-3.) */

/*       IWK = int work array of length at least 2*LWK. */

/*       LIST,LPTR,LEND = Data structure defining the trian- */
/*                        gulation.  Refer to Subroutine */
/*                        TRMESH. */

/* On output: */

/*       LWK = Number of arcs which intersect IN1-IN2 (but */
/*             not more than the input value of LWK) unless */
/*             IER = 1 or IER = 3.  LWK = 0 if and only if */
/*             IN1 and IN2 were adjacent (or LWK=0) on input. */

/*       IWK = Array containing the indexes of the endpoints */
/*             of the new arcs other than IN1-IN2 unless */
/*             IER > 0 or LWK = 0.  New arcs to the left of */
/*             IN1->IN2 are stored in the first K-1 columns */
/*             (left portion of IWK), column K contains */
/*             zeros, and new arcs to the right of IN1->IN2 */
/*             occupy columns K+1,...,LWK.  (K can be deter- */
/*             mined by searching IWK for the zeros.) */

/*       LIST,LPTR,LEND = Data structure updated if necessary */
/*                        to reflect the presence of an arc */
/*                        connecting IN1 and IN2 unless IER > */
/*                        0.  The data structure has been */
/*                        altered if IER >= 4. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if IN1 < 1, IN2 < 1, IN1 = IN2, */
/*                     or LWK < 0 on input. */
/*             IER = 2 if more space is required in IWK. */
/*                     Refer to LWK. */
/*             IER = 3 if IN1 and IN2 could not be connected */
/*                     due to either an invalid data struc- */
/*                     ture or collinear nodes (and floating */
/*                     point error). */
/*             IER = 4 if an error flag other than IER = 1 */
/*                     was returned by OPTIM. */
/*             IER = 5 if error flag 1 was returned by OPTIM. */
/*                     This is not necessarily an error, but */
/*                     the arcs other than IN1-IN2 may not */
/*                     be optimal. */

/*   An error message is written to the standard output unit */
/* in the case of IER = 3 or IER = 4. */

/* Modules required by EDGE:  LEFT, LSTPTR, OPTIM, SWAP, */
/*                              SWPTST */

/* Intrinsic function called by EDGE:  ABS */

/* *********************************************************** */


/* Local parameters: */

/* DPij =     Dot product <Ni,Nj> */
/* I =        DO-loop index and column index for IWK */
/* IERR =     Error flag returned by Subroutine OPTIM */
/* IWC =      IWK index between IWF and IWL -- NL->NR is */
/*              stored in IWK(1,IWC)->IWK(2,IWC) */
/* IWCP1 =    IWC + 1 */
/* IWEND =    Input or output value of LWK */
/* IWF =      IWK (column) index of the first (leftmost) arc */
/*              which intersects IN1->IN2 */
/* IWL =      IWK (column) index of the last (rightmost) are */
/*              which intersects IN1->IN2 */
/* LFT =      Flag used to determine if a swap results in the */
/*              new arc intersecting IN1-IN2 -- LFT = 0 iff */
/*              N0 = IN1, LFT = -1 implies N0 LEFT IN1->IN2, */
/*              and LFT = 1 implies N0 LEFT IN2->IN1 */
/* LP =       List pointer (index for LIST and LPTR) */
/* LP21 =     Unused parameter returned by SWAP */
/* LPL =      Pointer to the last neighbor of IN1 or NL */
/* N0 =       Neighbor of N1 or node opposite NR->NL */
/* N1,N2 =    Local copies of IN1 and IN2 */
/* N1FRST =   First neighbor of IN1 */
/* N1LST =    (Signed) last neighbor of IN1 */
/* NEXT =     Node opposite NL->NR */
/* NIT =      Flag or number of iterations employed by OPTIM */
/* NL,NR =    Endpoints of an arc which intersects IN1-IN2 */
/*              with NL LEFT IN1->IN2 */
/* X0,Y0,Z0 = Coordinates of N0 */
/* X1,Y1,Z1 = Coordinates of IN1 */
/* X2,Y2,Z2 = Coordinates of IN2 */


/* Store IN1, IN2, and LWK in local variables and test for */
/*   errors. */

    /* Parameter adjustments */
    --lend;
    --lptr;
    --list;
    iwk -= 3;
    --z__;
    --y;
    --x;

    /* Function Body */
    n1 = *in1;
    n2 = *in2;
    iwend = *lwk;
    if (n1 < 1 || n2 < 1 || n1 == n2 || iwend < 0) {
        goto L31;
    }

/* Test for N2 as a neighbor of N1.  LPL points to the last */
/*   neighbor of N1. */

    lpl = lend[n1];
    n0 = (i__1 = list[lpl], abs(i__1));
    lp = lpl;
L1:
    if (n0 == n2) {
        goto L30;
    }
    lp = lptr[lp];
    n0 = list[lp];
    if (lp != lpl) {
        goto L1;
    }

/* Initialize parameters. */

    iwl = 0;
    nit = 0;

/* Store the coordinates of N1 and N2. */

L2:
    x1 = x[n1];
    y1 = y[n1];
    z1 = z__[n1];
    x2 = x[n2];
    y2 = y[n2];
    z2 = z__[n2];

/* Set NR and NL to adjacent neighbors of N1 such that */
/*   NR LEFT N2->N1 and NL LEFT N1->N2, */
/*   (NR Forward N1->N2 or NL Forward N1->N2), and */
/*   (NR Forward N2->N1 or NL Forward N2->N1). */

/*   Initialization:  Set N1FRST and N1LST to the first and */
/*     (signed) last neighbors of N1, respectively, and */
/*     initialize NL to N1FRST. */

    lpl = lend[n1];
    n1lst = list[lpl];
    lp = lptr[lpl];
    n1frst = list[lp];
    nl = n1frst;
    if (n1lst < 0) {
        goto L4;
    }

/*   N1 is an interior node.  Set NL to the first candidate */
/*     for NR (NL LEFT N2->N1). */

L3:
    if (left_(&x2, &y2, &z2, &x1, &y1, &z1, &x[nl], &y[nl], &z__[nl])) {
        goto L4;
    }
    lp = lptr[lp];
    nl = list[lp];
    if (nl != n1frst) {
        goto L3;
    }

/*   All neighbors of N1 are strictly left of N1->N2. */

    goto L5;

/*   NL = LIST(LP) LEFT N2->N1.  Set NR to NL and NL to the */
/*     following neighbor of N1. */

L4:
    nr = nl;
    lp = lptr[lp];
    nl = (i__1 = list[lp], abs(i__1));
    if (left_(&x1, &y1, &z1, &x2, &y2, &z2, &x[nl], &y[nl], &z__[nl])) {

/*   NL LEFT N1->N2 and NR LEFT N2->N1.  The Forward tests */
/*     are employed to avoid an error associated with */
/*     collinear nodes. */

        dp12 = x1 * x2 + y1 * y2 + z1 * z2;
        dp1l = x1 * x[nl] + y1 * y[nl] + z1 * z__[nl];
        dp2l = x2 * x[nl] + y2 * y[nl] + z2 * z__[nl];
        dp1r = x1 * x[nr] + y1 * y[nr] + z1 * z__[nr];
        dp2r = x2 * x[nr] + y2 * y[nr] + z2 * z__[nr];
        if ((dp2l - dp12 * dp1l >= 0. || dp2r - dp12 * dp1r >= 0.) && (dp1l -
                dp12 * dp2l >= 0. || dp1r - dp12 * dp2r >= 0.)) {
            goto L6;
        }

/*   NL-NR does not intersect N1-N2.  However, there is */
/*     another candidate for the first arc if NL lies on */
/*     the line N1-N2. */

        if (! left_(&x2, &y2, &z2, &x1, &y1, &z1, &x[nl], &y[nl], &z__[nl])) {
            goto L5;
        }
    }

/*   Bottom of loop. */

    if (nl != n1frst) {
        goto L4;
    }

/* Either the triangulation is invalid or N1-N2 lies on the */
/*   convex hull boundary and an edge NR->NL (opposite N1 and */
/*   intersecting N1-N2) was not found due to floating point */
/*   error.  Try interchanging N1 and N2 -- NIT > 0 iff this */
/*   has already been done. */

L5:
    if (nit > 0) {
        goto L33;
    }
    nit = 1;
    n1 = n2;
    n2 = *in1;
    goto L2;

/* Store the ordered sequence of intersecting edges NL->NR in */
/*   IWK(1,IWL)->IWK(2,IWL). */

L6:
    ++iwl;
    if (iwl > iwend) {
        goto L32;
    }
    iwk[(iwl << 1) + 1] = nl;
    iwk[(iwl << 1) + 2] = nr;

/*   Set NEXT to the neighbor of NL which follows NR. */

    lpl = lend[nl];
    lp = lptr[lpl];

/*   Find NR as a neighbor of NL.  The search begins with */
/*     the first neighbor. */

L7:
    if (list[lp] == nr) {
        goto L8;
    }
    lp = lptr[lp];
    if (lp != lpl) {
        goto L7;
    }

/*   NR must be the last neighbor, and NL->NR cannot be a */
/*     boundary edge. */

    if (list[lp] != nr) {
        goto L33;
    }

/*   Set NEXT to the neighbor following NR, and test for */
/*     termination of the store loop. */

L8:
    lp = lptr[lp];
    next = (i__1 = list[lp], abs(i__1));
    if (next == n2) {
        goto L9;
    }

/*   Set NL or NR to NEXT. */

    if (left_(&x1, &y1, &z1, &x2, &y2, &z2, &x[next], &y[next], &z__[next])) {
        nl = next;
    } else {
        nr = next;
    }
    goto L6;

/* IWL is the number of arcs which intersect N1-N2. */
/*   Store LWK. */

L9:
    *lwk = iwl;
    iwend = iwl;

/* Initialize for edge swapping loop -- all possible swaps */
/*   are applied (even if the new arc again intersects */
/*   N1-N2), arcs to the left of N1->N2 are stored in the */
/*   left portion of IWK, and arcs to the right are stored in */
/*   the right portion.  IWF and IWL index the first and last */
/*   intersecting arcs. */

    iwf = 1;

/* Top of loop -- set N0 to N1 and NL->NR to the first edge. */
/*   IWC points to the arc currently being processed.  LFT */
/*   .LE. 0 iff N0 LEFT N1->N2. */

L10:
    lft = 0;
    n0 = n1;
    x0 = x1;
    y0 = y1;
    z0 = z1;
    nl = iwk[(iwf << 1) + 1];
    nr = iwk[(iwf << 1) + 2];
    iwc = iwf;

/*   Set NEXT to the node opposite NL->NR unless IWC is the */
/*     last arc. */

L11:
    if (iwc == iwl) {
        goto L21;
    }
    iwcp1 = iwc + 1;
    next = iwk[(iwcp1 << 1) + 1];
    if (next != nl) {
        goto L16;
    }
    next = iwk[(iwcp1 << 1) + 2];

/*   NEXT RIGHT N1->N2 and IWC .LT. IWL.  Test for a possible */
/*     swap. */

    if (! left_(&x0, &y0, &z0, &x[nr], &y[nr], &z__[nr], &x[next], &y[next], &
            z__[next])) {
        goto L14;
    }
    if (lft >= 0) {
        goto L12;
    }
    if (! left_(&x[nl], &y[nl], &z__[nl], &x0, &y0, &z0, &x[next], &y[next], &
            z__[next])) {
        goto L14;
    }

/*   Replace NL->NR with N0->NEXT. */

    swap_(&next, &n0, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    iwk[(iwc << 1) + 1] = n0;
    iwk[(iwc << 1) + 2] = next;
    goto L15;

/*   Swap NL-NR for N0-NEXT, shift columns IWC+1,...,IWL to */
/*     the left, and store N0-NEXT in the right portion of */
/*     IWK. */

L12:
    swap_(&next, &n0, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    i__1 = iwl;
    for (i__ = iwcp1; i__ <= i__1; ++i__) {
        iwk[(i__ - (1<<1)) + 1] = iwk[(i__ << 1) + 1];
        iwk[(i__ - (1<<1)) + 2] = iwk[(i__ << 1) + 2];
/* L13: */
    }
    iwk[(iwl << 1) + 1] = n0;
    iwk[(iwl << 1) + 2] = next;
    --iwl;
    nr = next;
    goto L11;

/*   A swap is not possible.  Set N0 to NR. */

L14:
    n0 = nr;
    x0 = x[n0];
    y0 = y[n0];
    z0 = z__[n0];
    lft = 1;

/*   Advance to the next arc. */

L15:
    nr = next;
    ++iwc;
    goto L11;

/*   NEXT LEFT N1->N2, NEXT .NE. N2, and IWC .LT. IWL. */
/*     Test for a possible swap. */

L16:
    if (! left_(&x[nl], &y[nl], &z__[nl], &x0, &y0, &z0, &x[next], &y[next], &
            z__[next])) {
        goto L19;
    }
    if (lft <= 0) {
        goto L17;
    }
    if (! left_(&x0, &y0, &z0, &x[nr], &y[nr], &z__[nr], &x[next], &y[next], &
            z__[next])) {
        goto L19;
    }

/*   Replace NL->NR with NEXT->N0. */

    swap_(&next, &n0, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    iwk[(iwc << 1) + 1] = next;
    iwk[(iwc << 1) + 2] = n0;
    goto L20;

/*   Swap NL-NR for N0-NEXT, shift columns IWF,...,IWC-1 to */
/*     the right, and store N0-NEXT in the left portion of */
/*     IWK. */

L17:
    swap_(&next, &n0, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    i__1 = iwf;
    for (i__ = iwc - 1; i__ >= i__1; --i__) {
        iwk[(i__ + (1<<1)) + 1] = iwk[(i__ << 1) + 1];
        iwk[(i__ + (1<<1)) + 2] = iwk[(i__ << 1) + 2];
/* L18: */
    }
    iwk[(iwf << 1) + 1] = n0;
    iwk[(iwf << 1) + 2] = next;
    ++iwf;
    goto L20;

/*   A swap is not possible.  Set N0 to NL. */

L19:
    n0 = nl;
    x0 = x[n0];
    y0 = y[n0];
    z0 = z__[n0];
    lft = -1;

/*   Advance to the next arc. */

L20:
    nl = next;
    ++iwc;
    goto L11;

/*   N2 is opposite NL->NR (IWC = IWL). */

L21:
    if (n0 == n1) {
        goto L24;
    }
    if (lft < 0) {
        goto L22;
    }

/*   N0 RIGHT N1->N2.  Test for a possible swap. */

    if (! left_(&x0, &y0, &z0, &x[nr], &y[nr], &z__[nr], &x2, &y2, &z2)) {
        goto L10;
    }

/*   Swap NL-NR for N0-N2 and store N0-N2 in the right */
/*     portion of IWK. */

    swap_(&n2, &n0, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    iwk[(iwl << 1) + 1] = n0;
    iwk[(iwl << 1) + 2] = n2;
    --iwl;
    goto L10;

/*   N0 LEFT N1->N2.  Test for a possible swap. */

L22:
    if (! left_(&x[nl], &y[nl], &z__[nl], &x0, &y0, &z0, &x2, &y2, &z2)) {
        goto L10;
    }

/*   Swap NL-NR for N0-N2, shift columns IWF,...,IWL-1 to the */
/*     right, and store N0-N2 in the left portion of IWK. */

    swap_(&n2, &n0, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    i__ = iwl;
L23:
    iwk[(i__ << 1) + 1] = iwk[(i__ - (1<<1)) + 1];
    iwk[(i__ << 1) + 2] = iwk[(i__ - (1<<1)) + 2];
    --i__;
    if (i__ > iwf) {
        goto L23;
    }
    iwk[(iwf << 1) + 1] = n0;
    iwk[(iwf << 1) + 2] = n2;
    ++iwf;
    goto L10;

/* IWF = IWC = IWL.  Swap out the last arc for N1-N2 and */
/*   store zeros in IWK. */

L24:
    swap_(&n2, &n1, &nl, &nr, &list[1], &lptr[1], &lend[1], &lp21);
    iwk[(iwc << 1) + 1] = 0;
    iwk[(iwc << 1) + 2] = 0;

/* Optimization procedure -- */

    *ier = 0;
    if (iwc > 1) {

/*   Optimize the set of new arcs to the left of IN1->IN2. */

        nit = iwc - (1<<2);
        i__1 = iwc - 1;
        optim_(&x[1], &y[1], &z__[1], &i__1, &list[1], &lptr[1], &lend[1], &
                nit, &iwk[3], &ierr);
        if (ierr != 0 && ierr != 1) {
            goto L34;
        }
        if (ierr == 1) {
            *ier = 5;
        }
    }
    if (iwc < iwend) {

/*   Optimize the set of new arcs to the right of IN1->IN2. */

        nit = iwend - (iwc<<2);
        i__1 = iwend - iwc;
        optim_(&x[1], &y[1], &z__[1], &i__1, &list[1], &lptr[1], &lend[1], &
                nit, &iwk[(iwc + (1<<1)) + 1], &ierr);
        if (ierr != 0 && ierr != 1) {
            goto L34;
        }
        if (ierr == 1) {
            goto L35;
        }
    }
    if (*ier == 5) {
        goto L35;
    }

/* Successful termination (IER = 0). */

    return 0;

/* IN1 and IN2 were adjacent on input. */

L30:
    *ier = 0;
    return 0;

/* Invalid input parameter. */

L31:
    *ier = 1;
    return 0;

/* Insufficient space reserved for IWK. */

L32:
    *ier = 2;
    return 0;

/* Invalid triangulation data structure or collinear nodes */
/*   on convex hull boundary. */

L33:
    *ier = 3;
/*      WRITE (*,130) IN1, IN2 */
/*  130 FORMAT (//5X,'*** Error in EDGE:  Invalid triangula', */
/*     .        'tion or null triangles on boundary'/ */
/*     .        9X,'IN1 =',I4,', IN2=',I4/) */
    return 0;

/* Error flag (other than 1) returned by OPTIM. */

L34:
    *ier = 4;
/*      WRITE (*,140) NIT, IERR */
/*  140 FORMAT (//5X,'*** Error in OPTIM (called from EDGE):', */
/*     .        '  NIT = ',I4,', IER = ',I1,' ***'/) */
    return 0;

/* Error flag 1 returned by OPTIM. */

L35:
    *ier = 5;
    return 0;
} /* edge_ */

int find_group	(	int	ix,
		int	iy,
		int	iz,
		int	grpid,
		EMData *	mg,
		EMData *	visited
	)

Definition at line 18696 of file util_sparx.cpp.

References EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), nx, ny, and EMAN::EMData::set_value_at().

Referenced by EMAN::Util::get_biggest_cluster().

{
        int offs[][3] = { {-1, 0, 0}, {1, 0, 0}, {0, -1, 0}, {0, 1, 0}, {0, 0, -1}, {0, 0, 1} };
        int noff = 6;

        int nx = visited->get_xsize();
        int ny = visited->get_ysize();
        int nz = visited->get_zsize();

        vector< point3d_t > pts;
        pts.push_back( point3d_t(ix, iy, iz) );
        visited->set_value_at( ix, iy, iz, (float)grpid );

        int start = 0;
        int end = pts.size();

        while( end > start ) {
                for(int i=start; i < end; ++i ) {
                        int ix = pts[i].x;
                        int iy = pts[i].y;
                        int iz = pts[i].z;

                        for( int j=0; j < noff; ++j ) {
                                int jx = ix + offs[j][0];
                                int jy = iy + offs[j][1];
                                int jz = iz + offs[j][2];

                                if( jx < 0 || jx >= nx ) continue;
                                if( jy < 0 || jy >= ny ) continue;
                                if( jz < 0 || jz >= nz ) continue;


                                if( (*mg)(jx, jy, jz)>0 && (*visited)(jx, jy, jz)==0.0 ) {
                                    pts.push_back( point3d_t(jx, jy, jz) );
                                    visited->set_value_at( jx, jy, jz, (float)grpid );
                                }

                        }
                }

                start = end;
                end = pts.size();
        }
        return pts.size();
}

int getnp_	(	double *	x,
		double *	y,
		double *	z__,
		int *	list,
		int *	lptr,
		int *	lend,
		int *	l,
		int *	npts,
		double *	df,
		int *	ier
	)

Definition at line 11872 of file util_sparx.cpp.

References abs.

{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    static int i__, n1;
    static double x1, y1, z1;
    static int nb, ni, lp, np, lm1;
    static double dnb, dnp;
    static int lpl;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/28/98 */

/*   Given a Delaunay triangulation of N nodes on the unit */
/* sphere and an array NPTS containing the indexes of L-1 */
/* nodes ordered by angular distance from NPTS(1), this sub- */
/* routine sets NPTS(L) to the index of the next node in the */
/* sequence -- the node, other than NPTS(1),...,NPTS(L-1), */
/* that is closest to NPTS(1).  Thus, the ordered sequence */
/* of K closest nodes to N1 (including N1) may be determined */
/* by K-1 calls to GETNP with NPTS(1) = N1 and L = 2,3,...,K */
/* for K .GE. 2. */

/*   The algorithm uses the property of a Delaunay triangula- */
/* tion that the K-th closest node to N1 is a neighbor of one */
/* of the K-1 closest nodes to N1. */


/* On input: */

/*       X,Y,Z = Arrays of length N containing the Cartesian */
/*               coordinates of the nodes. */

/*       LIST,LPTR,LEND = Triangulation data structure.  Re- */
/*                        fer to Subroutine TRMESH. */

/*       L = Number of nodes in the sequence on output.  2 */
/*           .LE. L .LE. N. */

/* The above parameters are not altered by this routine. */

/*       NPTS = Array of length .GE. L containing the indexes */
/*              of the L-1 closest nodes to NPTS(1) in the */
/*              first L-1 locations. */

/* On output: */

/*       NPTS = Array updated with the index of the L-th */
/*              closest node to NPTS(1) in position L unless */
/*              IER = 1. */

/*       DF = Value of an increasing function (negative cos- */
/*            ine) of the angular distance between NPTS(1) */
/*            and NPTS(L) unless IER = 1. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if L < 2. */

/* Modules required by GETNP:  None */

/* Intrinsic function called by GETNP:  ABS */

/* *********************************************************** */


/* Local parameters: */

/* DNB,DNP =  Negative cosines of the angular distances from */
/*              N1 to NB and to NP, respectively */
/* I =        NPTS index and DO-loop index */
/* LM1 =      L-1 */
/* LP =       LIST pointer of a neighbor of NI */
/* LPL =      Pointer to the last neighbor of NI */
/* N1 =       NPTS(1) */
/* NB =       Neighbor of NI and candidate for NP */
/* NI =       NPTS(I) */
/* NP =       Candidate for NPTS(L) */
/* X1,Y1,Z1 = Coordinates of N1 */

    /* Parameter adjustments */
    --x;
    --y;
    --z__;
    --list;
    --lptr;
    --lend;
    --npts;

    /* Function Body */
    lm1 = *l - 1;
    if (lm1 < 1) {
        goto L6;
    }
    *ier = 0;

/* Store N1 = NPTS(1) and mark the elements of NPTS. */

    n1 = npts[1];
    x1 = x[n1];
    y1 = y[n1];
    z1 = z__[n1];
    i__1 = lm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        ni = npts[i__];
        lend[ni] = -lend[ni];
/* L1: */
    }

/* Candidates for NP = NPTS(L) are the unmarked neighbors */
/*   of nodes in NPTS.  DNP is initially greater than -cos(PI) */
/*   (the maximum distance). */

    dnp = 2.;

/* Loop on nodes NI in NPTS. */

    i__1 = lm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        ni = npts[i__];
        lpl = -lend[ni];
        lp = lpl;

/* Loop on neighbors NB of NI. */

L2:
        nb = (i__2 = list[lp], abs(i__2));
        if (lend[nb] < 0) {
            goto L3;
        }

/* NB is an unmarked neighbor of NI.  Replace NP if NB is */
/*   closer to N1. */

        dnb = -(x[nb] * x1 + y[nb] * y1 + z__[nb] * z1);
        if (dnb >= dnp) {
            goto L3;
        }
        np = nb;
        dnp = dnb;
L3:
        lp = lptr[lp];
        if (lp != lpl) {
            goto L2;
        }
/* L4: */
    }
    npts[*l] = np;
    *df = dnp;

/* Unmark the elements of NPTS. */

    i__1 = lm1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        ni = npts[i__];
        lend[ni] = -lend[ni];
/* L5: */
    }
    return 0;

/* L is outside its valid range. */

L6:
    *ier = 1;
    return 0;
} /* getnp_ */

int i_dnnt

(

double *

x

)

Definition at line 7512 of file util_sparx.cpp.

Referenced by trplot_(), and vrplot_().

{
        return (int)(*x >= 0. ? floor(*x + .5) : -floor(.5 - *x));
}

int insert_	(	int *	k,
		int *	lp,
		int *	list,
		int *	lptr,
		int *	lnew
	)

Definition at line 12050 of file util_sparx.cpp.

Referenced by bdyadd_(), covsph_(), and intadd_().

{
    static int lsav;


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   07/17/96 */

/*   This subroutine inserts K as a neighbor of N1 following */
/* N2, where LP is the LIST pointer of N2 as a neighbor of */
/* N1.  Note that, if N2 is the last neighbor of N1, K will */
/* become the first neighbor (even if N1 is a boundary node). */

/*   This routine is identical to the similarly named routine */
/* in TRIPACK. */


/* On input: */

/*       K = Index of the node to be inserted. */

/*       LP = LIST pointer of N2 as a neighbor of N1. */

/* The above parameters are not altered by this routine. */

/*       LIST,LPTR,LNEW = Data structure defining the trian- */
/*                        gulation.  Refer to Subroutine */
/*                        TRMESH. */

/* On output: */

/*       LIST,LPTR,LNEW = Data structure updated with the */
/*                        addition of node K. */

/* Modules required by INSERT:  None */

/* *********************************************************** */


    /* Parameter adjustments */
    --lptr;
    --list;

    /* Function Body */
    lsav = lptr[*lp];
    lptr[*lp] = *lnew;
    list[*lnew] = *k;
    lptr[*lnew] = lsav;
    ++(*lnew);
    return 0;
} /* insert_ */

long int inside_	(	double *	p,
		int *	lv,
		double *	xv,
		double *	yv,
		double *	zv,
		int *	nv,
		int *	listv,
		int *	ier
	)

Definition at line 12109 of file util_sparx.cpp.

References b, ierr, intrsc_(), q, sqrt(), and TRUE_.

{
    /* Initialized data */

    static double eps = .001;

    /* System generated locals */
    int i__1;
    long int ret_val = 0;

    /* Builtin functions */
    //double sqrt(double);

    /* Local variables */
    static double b[3], d__;
    static int k, n;
    static double q[3];
    static int i1, i2, k0;
    static double v1[3], v2[3], cn[3], bp, bq;
    static int ni;
    static double pn[3], qn[3], vn[3];
    static int imx;
    static long int lft1, lft2, even;
    static int ierr;
    static long int pinr, qinr;
    static double qnrm, vnrm;
    extern /* Subroutine */ int intrsc_(double *, double *,
            double *, double *, int *);


/* *********************************************************** */

/*                                              From STRIPACK */
/*                                            Robert J. Renka */
/*                                  Dept. of Computer Science */
/*                                       Univ. of North Texas */
/*                                           renka@cs.unt.edu */
/*                                                   12/27/93 */

/*   This function locates a point P relative to a polygonal */
/* region R on the surface of the unit sphere, returning */
/* INSIDE = TRUE if and only if P is contained in R.  R is */
/* defined by a cyclically ordered sequence of vertices which */
/* form a positively-oriented simple closed curve.  Adjacent */
/* vertices need not be distinct but the curve must not be */
/* self-intersecting.  Also, while polygon edges are by defi- */
/* nition restricted to a single hemisphere, R is not so */
/* restricted.  Its interior is the region to the left as the */
/* vertices are traversed in order. */

/*   The algorithm consists of selecting a point Q in R and */
/* then finding all points at which the great circle defined */
/* by P and Q intersects the boundary of R.  P lies inside R */
/* if and only if there is an even number of intersection */
/* points between Q and P.  Q is taken to be a point immedi- */
/* ately to the left of a directed boundary edge -- the first */
/* one that results in no consistency-check failures. */

/*   If P is close to the polygon boundary, the problem is */
/* ill-conditioned and the decision may be incorrect.  Also, */
/* an incorrect decision may result from a poor choice of Q */
/* (if, for example, a boundary edge lies on the great cir- */
/* cle defined by P and Q).  A more reliable result could be */
/* obtained by a sequence of calls to INSIDE with the ver- */
/* tices cyclically permuted before each call (to alter the */
/* choice of Q). */


/* On input: */

/*       P = Array of length 3 containing the Cartesian */
/*           coordinates of the point (unit vector) to be */
/*           located. */

/*       LV = Length of arrays XV, YV, and ZV. */

/*       XV,YV,ZV = Arrays of length LV containing the Carte- */
/*                  sian coordinates of unit vectors (points */
/*                  on the unit sphere).  These values are */
/*                  not tested for validity. */

/*       NV = Number of vertices in the polygon.  3 .LE. NV */
/*            .LE. LV. */

/*       LISTV = Array of length NV containing the indexes */
/*               (for XV, YV, and ZV) of a cyclically-ordered */
/*               (and CCW-ordered) sequence of vertices that */
/*               define R.  The last vertex (indexed by */
/*               LISTV(NV)) is followed by the first (indexed */
/*               by LISTV(1)).  LISTV entries must be in the */
/*               range 1 to LV. */

/* Input parameters are not altered by this function. */

/* On output: */

/*       INSIDE = TRUE if and only if P lies inside R unless */
/*                IER .NE. 0, in which case the value is not */
/*                altered. */

/*       IER = Error indicator: */
/*             IER = 0 if no errors were encountered. */
/*             IER = 1 if LV or NV is outside its valid */
/*                     range. */
/*             IER = 2 if a LISTV entry is outside its valid */
/*                     range. */
/*             IER = 3 if the polygon boundary was found to */
/*                     be self-intersecti