simulation_4.mx3

Nx:= 1000
Ny:= 150 			
Nz:= 26
setGridSize(Nx, Ny, Nz)
setCellSize(4e-9, 4e-9, 4e-9)

Aex     = 13e-12
Xi      = 0.1
alpha   = 0.02 

distance:= 40
xoffset:= 20+37.5
nanoislandlength := 400

angle:= 360
ydim:=40e-9
zdim:=40e-9
translation := -(22.5 + 20)*1e-9
ytranslation:= (nanoislandlength + 40)*1e-9

bigR:= 20e-9			//radius of cylinder or our cuboid
currentdensity := 2*1e12
currentdensity2 := -2*1e12
currentdensity3 := -5.5*1e12

timestep := 100e-12
time1:= 8000e-12
time2:= 5000e-12
totaltime := 18000e-12

//Bxangle:= 1 //in degrees
//Bmag := 0.005 //Tesla

//current carrying wire
wireDW:=cuboid(4000e-9,ydim,zdim).transl(0,0,0)

//nanoisland circular ends
semicircle1:= cylinder(75e-9, 5e-9).transl(0, ((nanoislandlength-75)*0.5)*1e-9, 0)
semicircle2:= cylinder(75e-9, 5e-9).transl(0,-((nanoislandlength-75)*0.5)*1e-9, 0)
nanoisland:= cuboid(75e-9, (nanoislandlength-75)*1e-9, 5e-9).add(semicircle1).add(semicircle2)


//two nanoislands
wireSwitch1:=nanoisland.RotZ((360-angle)*pi/180).transl(300e-9,0,translation)
wireSwitch2:=nanoisland.RotZ((360-angle)*pi/180).transl(900e-9,0,translation)
wireSwitch3:=nanoisland.RotZ((360-angle)*pi/180).transl(-300e-9,0,translation)
wireSwitch4:=nanoisland.RotZ((360-angle)*pi/180).transl(-900e-9,0,translation)

wires:= wireSwitch1.add(wireSwitch2).add(wireSwitch3).add(wireSwitch4)

setGeom(wireDW.add(wires))
saveas(geom, "WireCross") 

//saveas(m, "geometry")

defRegion(1,wireDW)
defRegion(2, cuboid(40e-9,ydim,zdim).transl(-1980e-9,0,0))
defRegion(3, cuboid(40e-9,ydim,zdim).transl(1980e-9,0,0))
defRegion(4,wireSwitch1)
defRegion(5,wireSwitch2)
defRegion(6,wireSwitch3)
defRegion(7,wireSwitch4)

m.LoadFile("init_grid_cDW_2.ovf")

//m.setRegion(1, twodomain(1, 0, 0, 0, 0, -1, -1, 0, 0))
m.setRegion(2, uniform(1,0,0))
m.setRegion(3, uniform(-1,0,0))
m.setRegion(4, uniform(sin(angle*pi/180),cos(angle*pi/180), 0))
m.setRegion(5, uniform(sin(angle*pi/180),cos(angle*pi/180), 0))
m.setRegion(6, uniform(sin(angle*pi/180),cos(angle*pi/180), 0))
m.setRegion(7, uniform(sin(angle*pi/180),cos(angle*pi/180), 0))

Msat.setRegion(1, 800e3)
Msat.setRegion(2, 800e3)
Msat.setRegion(3, 800e3)
Msat.setRegion(4, 800e3)
Msat.setRegion(5, 800e3)
Msat.setRegion(6, 800e3)
Msat.setRegion(7, 800e3)

frozenspins.setRegion(2,1)
frozenspins.setRegion(3,1)

//saveas(m, "geometry2")
//saveas(B_ext, "beforeSC")

// Remove surface charges from left (mx=1) and right (mx=-1) sides to mimic infinitely long wire. We have to specify the region (0) at the boundaries.
BoundaryRegion := 1
MagLeft        := 1
MagRight       := -1
ext_rmSurfaceCharge(BoundaryRegion, MagLeft, MagRight)

saveas(m, "Init")
//saveas(B_ext, "afterSC")

//Minimize()

//saveas(m, "minimised")

// Schedule output
autosave(m, timestep)
tableadd(B_ext)
tableadd(B_ext.Region(1))
tableadd(B_ext.Region(4))
tableadd(B_ext.Region(5))
tableadd(B_ext.Region(6))
tableadd(B_ext.Region(7))
tableadd(m.Region(1))
tableadd(m.Region(4))
tableadd(m.Region(5))
tableadd(m.Region(6))
tableadd(m.Region(7))
tableadd(E_anis)
tableadd(E_demag)
tableadd(E_exch)
tableadd(E_therm)
tableadd(E_total)
tableadd(E_Zeeman)
tableautosave(timestep)

// Run the simulation with current through the sample

pol.setRegion(1, 0.56)
J.setRegion(1, vector(-1*currentdensity, 0, 0))


//Oersted field
//outside
maskout_1 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=0; j<70; j++{
		for k:=0; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_1.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_2 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=80; j<Ny; j++{
		for k:=0; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_2.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_3 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=18; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_3.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_4 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=0; k<8; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_4.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}





//inside
maskin := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=8; k<18; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) 		//radial distance of B field
			B := (mu0 * 1* currentdensity * ydim * zdim *d)/(2 * pi * bigR * bigR) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskin.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}




B_ext.add(maskout_1, 1)
B_ext.add(maskout_2, 1)
B_ext.add(maskout_3, 1)
B_ext.add(maskout_4, 1)
B_ext.add(maskin, 1)


saveas(B_ext, "oersted")


Run(time1)



//reset B_ext field
B_ext = vector (0,0,0)




// second current pulse
pol.setRegion(1, 0.56)
J.setRegion(1, vector(-1*currentdensity2, 0, 0))


//Oersted field
//outside
maskout_21 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=0; j<70; j++{
		for k:=0; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity2 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_21.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_22 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=80; j<Ny; j++{
		for k:=0; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity2 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_22.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_23 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=18; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity2 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_23.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_24 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=0; k<8; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity2 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_24.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}





//inside
maskin2 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=8; k<18; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) 		//radial distance of B field
			B := (mu0 * 1* currentdensity2 * ydim * zdim *d)/(2 * pi * bigR * bigR) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskin2.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}




B_ext.add(maskout_21, 1)
B_ext.add(maskout_22, 1)
B_ext.add(maskout_23, 1)
B_ext.add(maskout_24, 1)
B_ext.add(maskin2, 1)

saveas(B_ext, "oersted2")

Run(time2)


//reset B_ext field
B_ext = vector (0,0,0)


// second current pulse
pol.setRegion(1, 0.56)
J.setRegion(1, vector(-1*currentdensity3, 0, 0))


//Oersted field
//outside
maskout_31 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=0; j<70; j++{
		for k:=0; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity3 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_31.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_32 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=80; j<Ny; j++{
		for k:=0; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity3 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_32.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_33 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=18; k<Nz; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity3 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_33.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}

maskout_34 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=0; k<8; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) //radial distance of B field
			B := (mu0 * 1* currentdensity3 * ydim * zdim)/(2*pi*d) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskout_34.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}


//inside
maskin3 := newVectorMask(Nx, Ny, Nz)
for i:=0; i<Nx; i++{
	for j:=70; j<80; j++{
		for k:=8; k<18; k++{
			r := index2coord(i, j, k)
			x := r.X()
			y := r.Y()
			z := r.Z()
			d := sqrt(z*z + y*y) 		//radial distance of B field
			B := (mu0 * 1* currentdensity3 * ydim * zdim *d)/(2 * pi * bigR * bigR) 
			By := (-1*B*z)/d
			Bz := (B*y)/d
			maskin3.setVector(i, j, k, vector(0, By, Bz))
		}
	}
}


B_ext.add(maskout_31, 1)
B_ext.add(maskout_32, 1)
B_ext.add(maskout_33, 1)
B_ext.add(maskout_34, 1)
B_ext.add(maskin3, 1)


saveas(B_ext, "oersted3")


Run(totaltime - time1 - time2)