円を直線にする

  • こちらの続き
  • 円と直線にする
  • 曲率とその微分を図示しているわけだけれど、曲率とその微分とその積分と、曲線が閉じることについても考えたい
Cartesian2Spherical2<-function(x,y){
	if(x==0){list(x=x,y=0)}else{
		a<-y/x
		list(x=1/a*cos(4*atan(x)-pi/2),y=1/a*sin(4*atan(x)-pi/2)+1/a)
	}
}
Spherical2Cartesian2<-function(x,y){
	if(x^2+y^2==0){list(x=x,y=0)}else{
		a<-2*y/(x^2+y^2)
		if(a==0){list(x=x,y=y)}else{
			tmp<-atan(abs((y-1/a)/x))
			retx<-0
			rety<-0
			if(x>=0 & y <=1/a){
				tmp<-pi/2-tmp
			}else if(x>=0 & y> 1/a){
				tmp<-pi/2+tmp
			}else if (x<0 & y<=1/a){
				tmp<-3/2*pi+tmp
			}else{
				tmp<-3/2*pi-tmp
			}
			#tmpx<-x
			#tmpy<-y-1/a
			#tmp<--1/4*sign(a*x)*acos(a*x)+pi/2
			atmp<-tan(tmp/4)
			#print(a)
			list(x=atmp,y=a*atmp)
		}
	}
}

Nrep<-10
xssum<-NULL
col<-c()
Niter<-1000
ns<-sample(1:10,Nrep,replace=TRUE)
xmax<-5
for(rep in 1:Nrep){
	n<-ns[rep]
	col<-c(col,rep(rep,Niter))
	x<-y<-seq(from=0,to=xmax,length.out=Niter)
	as<-runif(n)
	#as<-c(1)
	for(i in 1:length(x)){
		y[i]<-sum(as*x[i]^(1:(n)))

	}

	x2<-x
	y2<-y
	for(i in 1:length(x2)){
		tmp<-Cartesian2Spherical2(x[i],y[i])
		x2[i]<-tmp$x
		y2[i]<-tmp$y
	}



	tmp<-matrix(c(x2,y2),ncol=2)
	xssum<-rbind(xssum,tmp)
}
for(i in 1:Nrep){
	plot(xssum[((i-1)*Niter+1):(i*Niter),1],xssum[((i-1)*Niter+1):(i*Niter),2],col=i,type="l",xlim=c(min(xssum),max(xssum)),ylim=c(min(xssum),max(xssum)))
	par(new=TRUE)

}
abline(0,1)

x2s<-xssum
for(i in 1:length(xssum[,1])){
	tmp<-Spherical2Cartesian2(xssum[i,1],xssum[i,2])
	x2s[i,]<-c(tmp$x,tmp$y)
}
for(i in 1:Nrep){
	plot(x2s[((i-1)*Niter+1):(i*Niter),1],x2s[((i-1)*Niter+1):(i*Niter),2],col=i,type="l",xlim=c(min(x2s[,1]),max(x2s[,1])),ylim=c(min(x2s[,2]),max(x2s[,2])))
	par(new=TRUE)

}
#abline(0,1)