LD := 70.0e-9
Lx:=500.0e-9
Ly:=Lx
Lz:=180e-9
nx:=256
ny:=256
nz:=64
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
dS:=3*dz
OuterCylinder:=Cylinder(Lx, Lz)
InnerCylinder:=Cylinder(Lx-2*dS, Lz-2*dS)
SkinLayer:=OuterCylinder.Sub(InnerCylinder)
defRegion(1, InnerCylinder)
defRegion(2, SkinLayer)
SetGeom(OuterCylinder)
EdgeSmooth=8
Dbulk.setRegion(1, 0.853e-3)
Dbulk.setRegion(2, 0.0e-9)


BD:= Dbulk.average()*Dbulk.average()/(2*Aex.average()*Msat.average())+mu0*Msat.average()
B:= 0.25
dB:=0.01
B_ext = vector(0, 0, B)


//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) + 2*PI + 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) 

m.setInShape(cylinder(3*LD  , Lz),uniform(0,0,1))
m.setInShape(cylinder(LD  , Lz).transl(LD,0,0),BlochSkyrmion(1, -1).transl(LD,0,0))
m.setInShape(cylinder(LD  , Lz).transl(-0.5*LD,-0.866*LD,0),BlochSkyrmion(1, -1).transl(-0.5*LD,-0.866*LD,0))
m.setInShape(cylinder(LD  , Lz).transl(-0.5*LD, 0.866*LD,0),BlochSkyrmion(1, -1).transl(-0.5*LD, 0.866*LD,0))
saveas(m,"ini.ovf")


MinimizerStop = 1e-5
minimize()
tablesave();
saveas(m,"ini250mT-mumax.ovf")


B=0.150
B_ext = vector(0, 0, B);
minimize();
saveas(m,"i150mT-mumax.ovf");

B=-0.050;
B_ext = vector(0, 0, B);
minimize();
tablesave();
saveas(m,"m050mT-mumax.ovf");

B=0.150;
B_ext = vector(0, 0, B);
minimize();
tablesave();
saveas(m,"p150mT-mumax.ovf");

B=0.200;
B_ext = vector(0, 0, B);
minimize();
tablesave();
saveas(m,"p200mT-mumax.ovf");

B=0.250;
B_ext = vector(0, 0, B);
minimize();
tablesave();
saveas(m,"p250mT-mumax.ovf");