Skip to content

Commit

Permalink
version 1.1-1
Browse files Browse the repository at this point in the history
  • Loading branch information
@gaborcsardi
Bin Wang authored and gaborcsardi committed Oct 31, 2011
1 parent 9e91363 commit b9fb23d
Show file tree
Hide file tree
Showing 7 changed files with 113 additions and 87 deletions.
6 changes: 3 additions & 3 deletions DESCRIPTION
View file
@@ -1,5 +1,5 @@
Package: spt
Version: 1.1
Version: 1.1-1
Date: 2011-10-31
Title: Sierpinski Pedal Triangle
Author: Bin Wang <bwang@jaguar1.usouthal.edu>.
Expand All @@ -8,6 +8,6 @@ Depends: R (>= 2.5.0), stats
Suggests: MASS
Description: This package collects algorithms related to Sierpinski pedal triangle (SPT).
License: Unlimited
Packaged: 2011-11-02 18:08:16 UTC; bwang
Packaged: 2011-11-05 12:59:00 UTC; bwang
Repository: CRAN
Date/Publication: 2011-11-02 19:03:20
Date/Publication: 2011-11-05 18:07:33
12 changes: 6 additions & 6 deletions MD5
View file
@@ -1,10 +1,10 @@
ca496720f39e42223e10686998a6b48c *DESCRIPTION
b50393939560bcc60c38ab94de897625 *DESCRIPTION
d008ac6798fdfeba4fb80de7176f68c3 *LICENCE
a9e19ab5a327f27fdd1ed031501c862f *NAMESPACE
5399d2f938beb63f730286f21262c7c8 *PORTING
bf3fb080f3e8ea2900c656f018266596 *R/kernel.R
2e46d215627025af7eed4bdd944b1b6e *NAMESPACE
180388a86e4e077d32202e9b1676be21 *PORTING
87bcb8da10ceb1042466d6fa5ecb2d12 *R/kernel.R
ee3784ad0fc5e1f68654fe3055a5481a *R/zzz.R
5f46927f84b4e710334ddf9942bafe34 *man/spt.Rd
4dbd0cb0e18afb5f8dc81a2cbf6b2b99 *man/spt.Rd
8290d2e9740414e315237f0d5d4024bb *src/Makevars
b0683d1aff464dec94cf1417d7e6871f *src/init.c
b6b2af4c6307892d3fdb0c5cfe0f6b66 *src/kernel.c
26da8718bb7ac51499b97b06da296311 *src/kernel.c
2 changes: 1 addition & 1 deletion NAMESPACE
View file
@@ -1,4 +1,4 @@
useDynLib(spt, .registration = TRUE, .fixes = "F_")
useDynLib(spt, .registration = TRUE, .fixes = ".F_")
export(.sptConnect)
##export everything that does not start with a .
exportPattern("^[^\\.]")
Expand Down
6 changes: 5 additions & 1 deletion PORTING
View file
Expand Up @@ -10,4 +10,8 @@ option chaosgame=FALSE to determine how to draw the SPT. 4) iteration
number will not be specified. Taking default value 10000 for
"chaosgame" and 15 for regular drawing. 5) removing "main=NULL". It
can be specified, but not prespecified. 6) "tol=0.00001" is removed,
and taken as default value.
and taken as default value.

11/4/11: In this version, (1) we fixed a bug in computing the
SPT dimention. (2) We restricted the expiration of the external C/Fortran
functions to the namespace. (3) output a class 'sped' for further study.
131 changes: 76 additions & 55 deletions R/kernel.R
View file
@@ -1,10 +1,15 @@
spt <- function(A,B,x0,y0,iter,
spt <- function(A,B,x0,y0,iter,plot=TRUE,
spt=TRUE,chaosgame=FALSE,main=NULL,tol=0.0001,...){
if(missing(A)||missing(B))
stop("Sizes of at least two angles need to be specified.")
if(A<=0 || B<=0 || A>=180 || B >=180 || A+B>=180)
stop("Invaid angle size(s) specified: 0<A<180, 0<B<180, and 0<A+B<180.")
A1 = A; B1 = B; C1 = 180-A-B
A1 = A; B1 = B; C1 = 180.0-(A1+B1)
triname = ifelse(spt,"SPT","ST")
triname = paste(triname,"(",as.character(format(A1,20)),",",
as.character(format(B1,20)),",",
as.character(format(180-A1-B1)),")",sep='')

angles = sort(c(A1,B1,C1),decreasing=TRUE)
minA = min(angles); maxA = max(angles);
r = pi/180; h=100;
Expand All @@ -14,66 +19,82 @@ spt <- function(A,B,x0,y0,iter,
ABC = list(A=A,B=B,C=C,xa=xa,xb=xb,xc=xc,ya=ya,yb=yb,yc=yc,A0=A0,B0=B0,C0=C0);

xmin = min(xa,xb,xc); xmax=max(xa,xb,xc);
ymin = 0; ymax=h;
triname = ifelse(spt,"spt","pt")
sptdim = ifelse(maxA>90,NA,
round(.Fortran(F_SptDimFP, as.double(angles),as.double(0))[[2]],3))
if(is.null(main)){
main=paste(triname,"(",A1,",",B1,",",180-A1-B1,")",sep='')
if(is.na(sptdim))
main=paste(main, ", Dim = ",sptdim,sep='')
}
if(maxA > 90 && spt){
delta = abs(h * sin(maxA))*.25
xmin = xmin - delta;
delta = (xmax-xmin-ymax+ymin)*.75;
ymin = ymin - delta;
ymax = ymax + delta;
}

## draw the initial triangle
plot(0,0, type='n', bty='n', xaxt='n',yaxt='n',
xlab='', ylab='', asp=1, main = main,
xlim=c(xmin,xmax), ylim=c(ymin,ymax))
lines(c(ABC$xa,ABC$xb,ABC$xc,ABC$xa),
c(ABC$ya,ABC$yb,ABC$yc,ABC$ya),
lwd=2);

if(chaosgame){
if(missing(x0)) x0 = (xmax-xmin)*runif(1)+xmin;
if(missing(y0)) y0 = (ymax-ymin)*runif(1)+ymin;
if(missing(iter)) iter=10000;

points(x0,y0,pch='*', col=5)
points(ABC$xa,ABC$ya,pch=20, col=2)
points(ABC$xb,ABC$yb,pch=20, col=3)
points(ABC$xc,ABC$yc,pch=20, col=4)
if(spt){
.GameSPT(x0,y0,ABC,iter)
}else{
.GameST(x0,y0,ABC,iter)
ymin = 0; ymax=h;
if(maxA>90) sptdim = NA
else if(maxA==90) sptdim = 2.0
else sptdim = .Fortran(.F_SptDimFP, as.double(angles),as.double(0))[[2]]
## sptdim = .Fortran(F_SptDimFP, as.double(angles),as.double(0))[[2]]
if(plot){
if(is.null(main)){
## if(is.na(sptdim)||sptdim < 0)
cond1 = is.null(sptdim)
cond2 = is.na(sptdim)
if(cond1 || cond2)
main=paste(triname, ", Dim = NA",sep='')
else
main=paste(triname, ", Dim = ",as.character(format(sptdim,20)),sep='')
}
}else{
if(missing(iter)) iter = 20
Tri0 = cbind(c(xa,xb,xc), c(ya,yb,yc));
tol = tol * (xmax-xmin) * (ymax-ymin)
if(iter>0){
if(!spt) .ptpedal(iter,Tri0,tol)
else if(angles[1]==90){
.sptpedal2(iter,Tri0,tol)
if(maxA > 90 && spt){
delta = abs(h * sin(maxA))*.25
xmin = xmin - delta;
delta = (xmax-xmin-ymax+ymin)*.75;
ymin = ymin - delta;
ymax = ymax + delta;
}

## draw the initial triangle
plot(0,0, type='n', bty='n', xaxt='n',yaxt='n',
xlab='', ylab='', asp=1, main = main,
xlim=c(xmin,xmax), ylim=c(ymin,ymax))
lines(c(ABC$xa,ABC$xb,ABC$xc,ABC$xa),
c(ABC$ya,ABC$yb,ABC$yc,ABC$ya),
lwd=2);

if(chaosgame){
if(missing(x0)) x0 = (xmax-xmin)*runif(1)+xmin;
if(missing(y0)) y0 = (ymax-ymin)*runif(1)+ymin;
if(missing(iter)) iter=10000;

points(x0,y0,pch='*', col=5)
points(ABC$xa,ABC$ya,pch=20, col=2)
points(ABC$xb,ABC$yb,pch=20, col=3)
points(ABC$xc,ABC$yc,pch=20, col=4)
if(spt){
.GameSPT(x0,y0,ABC,iter)
}else{
iter = ifelse(maxA>90, min(12,iter),iter)
.sptpedal3(iter,Tri0,tol)
.GameST(x0,y0,ABC,iter)
}
}else{
if(missing(iter)) iter = 20
Tri0 = cbind(c(xa,xb,xc), c(ya,yb,yc));
tol = tol * (xmax-xmin) * (ymax-ymin)
if(iter>0){
if(!spt) .ptpedal(iter,Tri0,tol)
else if(angles[1]==90){
.sptpedal2(iter,Tri0,tol)
}else{
iter = ifelse(maxA>90, min(12,iter),iter)
.sptpedal3(iter,Tri0,tol)
}
}
}
}
invisible(round(sptdim,3));

if(spt){
out = structure(list(angles,
data.name = triname,
Dim = sptdim, iter=iter
), class = "spt")
cat("\n\nThe SPT dimension is", format(sptdim,10),"\n\n");
}else out=NULL

invisible(out)
}

## Draw PT/SPT: tri1 is the initial triangle, n is the iteration number.
.ptpedal <- function(n, tri1,tol){
if(n>1){
tmp = .Fortran(F_PTChild, as.double(as.vector(tri1)), as.double(rep(0,6)));
tmp = .Fortran(.F_PTChild, as.double(as.vector(tri1)), as.double(rep(0,6)));
result = matrix(tmp[[2]], ncol=2, nrow=3);
polygon(result[,1],result[,2])
tri1.1 = rbind(tri1[3,], result[2,], result[1,]);
Expand All @@ -89,7 +110,7 @@ spt <- function(A,B,x0,y0,iter,

.sptpedal2 <- function(n, tri1,tol){
if(n>1){
tmp = .Fortran(F_SPTChild2, as.double(as.vector(tri1)), as.double(rep(0,2)));
tmp = .Fortran(.F_SPTChild2, as.double(as.vector(tri1)), as.double(rep(0,2)));
tmp = matrix(tmp[[2]],nrow=1);
segments(tmp[1],tmp[2],tri1[1,1],tri1[1,2]);
tri1.1 = rbind(tmp, tri1[3,], tri1[1,]);
Expand All @@ -103,7 +124,7 @@ spt <- function(A,B,x0,y0,iter,

.sptpedal3 <- function(n, tri1, tol){
if(n>1){
tmp = .Fortran(F_SPTChild3, as.double(as.vector(tri1)), as.double(rep(0,6)));
tmp = .Fortran(.F_SPTChild3, as.double(as.vector(tri1)), as.double(rep(0,6)));
result = matrix(tmp[[2]], ncol=2, nrow=3);
polygon(result[,1],result[,2])
tri1.1 = rbind(tri1[3,], result[2,], result[1,]);
Expand All @@ -119,7 +140,7 @@ spt <- function(A,B,x0,y0,iter,

.Stop <- function(triangle){
axis = as.vector(triangle);
.Fortran(F_stri, as.double(axis), as.double(0))[[2]]
.Fortran(.F_stri, as.double(axis), as.double(0))[[2]]
}

##########################################################################
Expand Down
12 changes: 4 additions & 8 deletions man/spt.Rd
View file
Expand Up @@ -3,23 +3,19 @@

\name{spt}
\alias{spt}
\alias{F_PTChild}
\alias{F_SPTChild2}
\alias{F_SPTChild3}
\alias{F_SptDimFP}
\alias{F_stri}

\title{Sierpinski Pedal Triangle and Sierpinski Triangle}
\description{
To draw Sierpinski pedal triangles (SPTs) and Sierpinski triangles (PTs).
}
\usage{
spt(A,B,x0,y0,iter, spt=TRUE,chaosgame=FALSE,main=NULL,tol=0.0001,...)
spt(A,B,x0,y0,iter,plot=TRUE, spt=TRUE,chaosgame=FALSE,main=NULL,tol=0.0001,...)
}
\arguments{
\item{A,B}{The degrees of two of the three angles of a triangle.}
\item{x0,y0}{The initial position of the point if "chaosgame=TRUE".}
\item{iter}{Iteration number in drawing the SPTs/PTs.}
\item{plot}{Draw the SPT/PT (default) or compute the SPT dimention alone.}
\item{spt}{SPT ("spt=TRUE") or PT ("spt=FALSE").}
\item{chaosgame}{To draw the SPTs/PTs using regular method or via a chaos game. }
\item{main}{Specify the title of the plot.}
Expand Down Expand Up @@ -53,8 +49,8 @@

spt(45,55,spt=FALSE)
spt(45,55,spt=TRUE)
spt(120,30,iter=15)
spt(10,10,iter=10)
spt(120,30,iter=3)
spt(20,20,iter=7)

}
\keyword{stats}
Expand Down
31 changes: 18 additions & 13 deletions src/kernel.c
View file
Expand Up @@ -89,7 +89,7 @@ void PTChild(double *xy, double *result)
}

double fsptdim(double a, double b, double c, double d){
return(pow(cos(a), d)+pow(cos(b),d)+pow(cos(c),d)-1.);
return(R_pow(cos(a), d)+R_pow(cos(b),d)+R_pow(cos(c),d)-1.);
}

void SptDimFP(double *angles, double *dim){
Expand All @@ -101,21 +101,26 @@ void SptDimFP(double *angles, double *dim){
Secant method is used. (a generalization of Secant method:
Method of False position)
*/
double a,b,c;
double x0=1., x1 = 3., xi, fa,fb,fc = 1.; //initial values;
double a,b,c, tol = 6.123234e-17;
double x0=log(3.0)/log(2.0), x1 = 2.0, xi, fa,fb,fc = 1.; //initial values;
a = angles[0]*PI/180.;
b = angles[1]*PI/180.;
c = angles[2]*PI/180.;
fa = fsptdim(a,b,c,x0); fb= fsptdim(a,b,c,x1);
while(fabs(fc) >.0000000001 && fa * fb <0.){
xi = (x1 * fa - x0 * fb)/(fa-fb);
fc = fsptdim(a,b,c,xi);
if(fa*fc<0){
x1=xi;}
else{
x0=xi;
}
//a suppose to be the larget angle
if(R_pow(cos(a),x0)<= tol)
dim[0] = 2.0;
else{
fa = fsptdim(a,b,c,x0); fb= fsptdim(a,b,c,x1);
while(fabs(fc) > tol){
xi = 0.5 *(x0+x1);
fc = fsptdim(a,b,c,xi);
if(fc<0){
x1=xi;}
else{
x0=xi;
}
fa = fsptdim(a,b,c,x0); fb= fsptdim(a,b,c,x1);
}
dim[0] = xi;
}
dim[0] = xi;
}

0 comments on commit b9fb23d

Please sign in to comment.