#Numerical exercise for vector algebra, vorticity and Gauss and Stokes law

#This is just a definition of a function to plot vectorplots, you do not have to understand it...
par.uin<-function() 
  # determine scale of inches/userunits in x and y
  # from http://tolstoy.newcastle.edu.au/R/help/01c/2714.html
  # Brian Ripley Tue 20 Nov 2001 - 20:13:52 EST
 {
    u <- par("usr") 
    p <- par("pin") 
    c(p[1]/(u[2] - u[1]), p[2]/(u[4] - u[3]))
 }

quiver<-function(lon,lat,u,v,scale=1,length=0.2,maxv=max(abs(na.omit(u)),abs(na.omit(v))), ...) 
 # first stab at matlab's quiver in R
 # from http://tolstoy.newcastle.edu.au/R/help/01c/2711.html
 # Robin Hankin Tue 20 Nov 2001 - 13:10:28 EST
  {
    ypos <- lat[col(u)] 
    xpos <- lon[row(u)]

    speed <- sqrt(u*u+v*v) 

    u <- u*scale/maxv
    v <- v*scale/maxv
    

    matplot(xpos,ypos,type="p",cex=0, ...) 
    arrows(xpos,ypos,xpos+u,ypos+v,length=length*min(par.uin())) 
  }


#Differentiation of a scalar 2D field s[x,y] in x direction
#It assumes that the distance between the gridpoints in each direction is 1 (h=1)
diff.x<-function(s)
{
	Nx<-dim(s)[1]   #Size of the field in X direction
	Ny<-dim(s)[2]   #Size of the field in Y direction

	h<-1            #h = 1
	
	result<-matrix(NA,Nx,Ny) #Definiton a Nx x Ny field, filled with NA's = missing values
	
	for (ix in 2:(Nx-1))   #loop through all cells 
		for (iy in 1:Ny)
		{
			result[ix,iy]<-(s[ix+1,iy]-s[ix-1,iy])/(2*h)  #Calculate the symetric derivative
		}

	return(result)
}


#Differentiation of a scalar 2D field s[x,y] in y direction
#It assumes that the distance between the gridpoints in each direction is 1 (h=1)
diff.y<-function(s)
{
	Nx<-dim(s)[1]   #Size of the field in X direction
	Ny<-dim(s)[2]   #Size of the field in Y direction

	h<-1            #h = 1
	
	result<-matrix(NA,Nx,Ny) #Definiton a Nx x Ny field, filled with NA's = missing values
	
	for (ix in 2:Nx)   #loop through all cells 
		for (iy in 2:(Ny-1))
		{
			result[ix,iy]<-(s[ix,iy+1]-s[ix,iy-1])/(2*h)  #Calculate the symetric derivative
		}

	return(result)
}



#2D Area Integral = line integral of the normal component
#around a rectangular box with the index given by the index of the lower left corner
#i.x1,i.y1 and the upper right corner i.x2,i.y2
#Again it is assumed that the grid distance = 1

AreaIntegral<-function(vx,vy,i.x1,i.y1,i.x2,i.y2)
{

#sum up the normal component: right on top side + down on right side + left on lower side + up on the left side
	integral<-sum(vy[i.x1:i.x2,i.y2])+sum(vx[i.x2,i.y2:i.y1])-sum(vy[i.x2:i.x1,i.y1])-sum(vx[i.x1,i.y1:i.y2])       
	
	return(integral)
}

#Line Integral around a rectangular box with the index given by the index of the lower left corner
#i.x1,i.y1 and the upper right corner i.x2,i.y2, direction is clockwise


LineIntegral<-function(vx,vy,i.x1,i.y1,i.x2,i.y2)
{
#sum up the tangential component: right on top side + down on right side + left on lower side + up on the left side
	integral<-sum(vx[i.x1:i.x2,i.y2])-sum(vy[i.x2,i.y2:i.y1])-sum(vx[i.x2:i.x1,i.y1])+sum(vy[i.x1,i.y1:i.y2])       
	return(integral)
}




#2D Volume integral of a rectangular box with the index given by the index of the lower left corner
#i.x1,i.y1 and the upper right corner i.x2,i.y2

VolumeIntegral<-function(s,i.x1,i.y1,i.x2,i.y2)
{
	integral<-sum(s[i.x1:i.x2,i.y1:i.y2])
	return(integral)
}


#Conversion from Cartesian to Polar coordinates
#This uses atan2: a two-argument function that computes 
# the arctangent of y / x given y and x, but with a range of ( - pi,pi]. 
#for more information, e.g. http://en.wikipedia.org/wiki/Atan2

polar<-function(x,y)
{
	r<-sqrt(x^2+y^2)
	phi<-atan2(y,x)
	return(c(r,phi))
}


#Conversion from Polar to Cartesian coordinates

cartesian<-function(r,phi)
{
	x<-r*cos(phi)
	y<-r*sin(phi)
	return(c(x,y))
}

#####################################
############################### Start of main program #####################


x<-(-10):10   #Define the midpoints of the gridboxes
y<-(-10):10

vx<-matrix(0,21,21)  #Define fields for the x and y components 
vy<-matrix(0,21,21)  #of the vector field


#Initalize a rotating disc
omega=0.1 #angular speed

for (iX in 1:21) for (iY in 1:21)
{
	
	pos.polar<-polar(x[iX],y[iY]) #Get the actual position in polar coordinates

	#Calculate the velocity vector in Cartesian coordinates which is just the difference between
      # (r,phi) and (r,phi+omega)
	t<-cartesian(pos.polar[1],pos.polar[2])-cartesian(pos.polar[1],pos.polar[2]+omega)

	#put the result in the velocity vectors
	vx[iX,iY]<-t[1]	
	vy[iX,iY]<-t[2]
}

quiver(x,y,vx,vy)