Lx := 1000.0e-9
Ly := 1000.0e-9
Lz := 180.0e-9
nx := 400
ny := 400
nz := 72
dx := Lx/nx
dy := Ly/ny
dz := Lz/nz
SetMesh(nx,ny,nz,dx,dy,dz,0,0,0)

Msat = 384e3
Aex = 4.75e-12
Dbulk = 0.853e-3//4*pi*Aex.average()/LD
B:= 0.243
B_ext = vector(0, 0, B)
BD:= Dbulk.average()*Dbulk.average()/(2*Aex.average()*Msat.average())+mu0*Msat.average()
LD := 70.0e-9//4*pi*Aex.average()/Dbulk.average() 
//Setting the geometry
dS:=3*dx//thickness of the FIB damaged layer
OuterBox:=Cuboid(Lx, Ly, Lz)
InBox:=Cuboid(Lx-dS, Ly-dS, Lz-dS)
SkinLayer:=OuterBox.Sub(InBox)
defRegion(1, InBox)
defRegion(2, SkinLayer)
Dbulk.setRegion(1, Dbulk.average())
Dbulk.setRegion(2, 0.0)



TableAdd(Edens_total)

//Setting initial state
//setting vortex-cone 
vortex_cone := newVectorMask(nx, ny, nz)
for k:=0; k<Nz; k++{
    for j:=0; j<Ny; j++{
        for i:=0; i<Nx; i++{
            r := index2coord(i, j, k)
            x := r.X()
            y := r.Y()
            z := r.Z()
            T := acos(B/BD)
            F := atan2(y,x) + PI/2.0 + 2*PI*z/LD
            vortex_cone.setVector(i, j, k, vector(sin(T)*cos(F), sin(T)*sin(F), cos(T)))
        }
    }
    print("setting layer:", k," of ", nz)
}
m.SetArray(vortex_cone) 
saveas(m,"vortex_cone_0p5pi.ovf")


// Since the methods SetVector and SetCell are extrimly slow
// it is more efficient to save the vortex-cone state in ovf file
// m.LoadFile("vortex_cone.ovf")

// Setting initial state, see Extended Data Fig. 8
d1 := 4*LD
d2 := 3*LD
d3 := 2*LD
h1 := Lz/3.5 //state 1
// h1 := Lz //state 2
// the outer ring:
m.SetInShape(Cylinder(d1, Lz), uniform(0,0,1))//magnetization up
m.SetInShape(Cylinder(d2, h1).Sub(Cylinder(d3, h1)), uniform(0,0,-1))//magnetization down
// the skyrmion string 
m.SetInShape(Cylinder(LD, Lz), uniform(0,0,-1))//magnetization up
//m.SetInShape(Cylinder(LD, Lz), BlochSkyrmion(1, -1))//optional 
saveas(m,"ini.ovf")

// relax()//might be required for some initial configurations 
minimize()
save(m)