f<-function(w,d)

{

  l<-w/sqrt(w^2+1)*exp(-atan(w)/w)*exp(-d/w)

  r<-sin(d)

  return(l-r)

}

w<-seq(0.01,60,0.01)

d<-rep(0,length(w))

s<-rep(0,length(w))

for(i in 1:length(w))

{

  root<-uniroot(f,c(0,pi/2),w=w[i],tol=0.01)

  d[i]<-root$root

  s[i]<-pi-d[i]-atan(w[i])

}

library(ggplot2)

DataSet1<-data.frame(w,value=d,lab=rep("d",length(w)))

DataSet2<-data.frame(w,value=s,lab=rep("s",length(w)))

DataSet <- rbind(DataSet1,DataSet2)

p<-ggplot(data=DataSet,aes(w,value,color=lab))

p+geom_line()+

  scale_colour_hue("value",breaks=c("d","s"),labels=c(expression(delta),expression(theta)))+

  labs(x=(expression(omega*RC/rad)),y=expression(list(delta,theta)/rad))+

  xlim(0,60)+ylim(0,pi)+

  scale_x_continuous(expand = c(0,0))+

  scale_y_continuous(breaks=round(c(0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6,pi),digits=2),

                     labels=expression(0,pi/6,pi/3,pi/2,2*pi/3,5*pi/6,pi),expand = c(0, 0))

AngleToRadian<-function(alpha)

{

  return(alpha*pi/180)

}

RadianToAngle<-function(alpha)

{

  return(alpha/pi*180)

}

f<-function(alpha,theta,varphi)

{

  var1<-sin(alpha+theta-varphi)

  if(varphi==0)

    var2<-0

  else

    var2<-sin(alpha-varphi)*exp(-theta/tan(varphi))

  return(var1-var2)

}

alpha<-seq(0,180-2,0.5)

varphi<-seq(0,90,15)

data<-data.frame()

for(a in alpha)

{

  for(v in varphi)

  {

    if(v<=a){

    t=uniroot(f,c(1e-3,pi+1e-3),alpha=AngleToRadian(a),

              varphi=AngleToRadian(v),tol=0.001)

    result<-c(a,RadianToAngle(t$root),v)

    data<-rbind(data,result)}

  }

}

names(data)<-c("alpha","theta","varphi")

data$I_VTN<-rep(0,dim(data)[1])

for(i in 1:dim(data)[1])

{

  a<-AngleToRadian(data$alpha[i])

  t<-AngleToRadian(data$theta[i])

  if(data$varphi[i]==90)

    v<-AngleToRadian(data$varphi[i]-1)

  else

    v<-AngleToRadian(data$varphi[i])

  data$I_VTN[i]<-1/(2*sqrt(pi))*sqrt(t-sin(t)*cos(2*a+v+t)/cos(v))

}

library(ggplot2)

theme_set(theme_light())

ggplot(data=data,aes(alpha,theta,color=factor(varphi)))+

  geom_line()+scale_colour_hue(expression(varphi/(degree)))+

  labs(x=(expression(alpha/(degree))),y=expression(theta/(degree)))

ggplot(data=data,aes(alpha,I_VTN,color=factor(varphi)))+geom_line()+

  geom_line()+scale_colour_hue(expression(varphi/(degree)))+

  labs(x=(expression(alpha/(degree))),y=expression(I[VTN]))