//Make a Sky in a SAF simple

SnapshotFormat = "png";
outputformat = OVF1_TEXT

// Set the Mesh

repeat_units:=2 //number of magnetic layers

cellsizeXY := 1.5625e-9
cellsizeZ := 0.4e-9 
Nx:=256; Ny:=Nx 
SetPBC(5,5,0)
SetGridSize(Nx,Ny,2*repeat_units)
SetCellSize(cellsizeXY, cellsizeXY, cellsizeZ)
// for the mesh we have 2X the magnetic layer -> for a simple SAF : 4
//I think for a simple SAF, it is possible to make with just 3 layer mesh (but I didn't run that code in thoses conditions)

// Set up material parameters (& external fields if u want)



Msat_value:= 580e3 //It is important to define that quantity as a simple Scalair for the following
//If you whant to un-balance the SAF you need also :
Ms1i := Msat_value
Ms2i := Msat_value


Msat = Msat_value
Aex=15e-12
Dind=3e-3
Ku1 = 0.7e6
alpha = 0.01
AnisU = vector(0,0,1)
B_ext = vector(0, 0, 0)


//define the geometry (1/2) layer of the mesh are magnetic
//define the regions coresponding to each magnetic layer

MagLayers:= Layer(0) // Initialise shape
for i:=0 ; i < repeat_units ; i++ {
    MagLayers=MagLayers.add(Layer(2*i)) // Select every second layer
	DefRegion(i, Layer(2*i)) //for later references to layers
}
SetGeom(MagLayers) // Make these layers magnetic

//////////////////////
// un-balenced SAF ///
//////////////////////

for i:=0 ; i < repeat_units ; i+=2 {
Msat.SetRegion(i, Ms1i)
Msat.SetRegion(1+i, Ms2i)
}

////////////////////////////
////Set up the RKKY////
////////////////////////////
RKKYAFM:= -0.00005

B_ext = vector(0,0,0)
TableAdd(B_ext)
TableAdd(m.region(0))
TableAdd(m.region(1))
TableAddVar(RKKYAFM, "RKKY", "j/m^2")

////////////////////////////
//Make the RKKY as a field//
////////////////////////////


for i:=0 ; i <= repeat_units-2 ; i++ {
	B_AFM1 := Masked(Mul(Const(RKKYAFM/(cellsizeZ*Ms1i)),shifted(m,0,0,-2)),layer(2*i))
	B_AFM2 := Masked(Mul(Const(RKKYAFM/(cellsizeZ*Ms2i)),shifted(m,0,0,2)),layer(2*i+2))
	B_AFM := Add(B_AFM1,B_AFM2)
	Edens_AFM := Mul(Msat,Dot(Mul(Const(-1),m),B_AFM))
	AddFieldTerm(B_AFM)
	AddEdensTerm(Edens_AFM)
}



////////////////
//Show-time !///
////////////////



for i:=0 ; i < repeat_units ; i++ {
		m.SetRegion(i, NeelSkyrmion(pow(-1,i), pow(-1,i+1)).scale(3,3,1))
		//m.setRegion(i, randomMag())
	}

saveas(m, "initialised_m")	
relax()
minimize()

saveas(m, "initial_m")

H:=0.5
B_ext = vector(H,0,50e-3)
relax()

A:=20e-3
f:=88.0e9


B_ext = vector(H,0,50e-3 + A*sin(2*pi*f*t))
tableautosave(2e-12)
Autosave(m,0.5e-10)
run(10e-9)
