GOSUALISATION.go

package main

import (
	"fmt"
	"math"
	"myfunc"
	"os"
	"strconv"
)

//Фрагмент программы, имитирующий omf-файл. Исходная намагниченность
//рассчитана на сетке 128x32x256, массивы fila, turn, turn имеют
//размеры 126х30х254, после "прореживания" пишем в omf-файл массив
//с размерами 63х15х127

func main() {

	var Mtrxd2 = map[int]map[int]float64{}
	var valx = myfunc.Range(1, 256)
	var valy = myfunc.Range(1, 256)
	for _, i := range valx {
		Mtrxd2[i] = make(map[int]float64)
		for _, j := range valy {
			Mtrxd2[i][j] = 0
		}
	}
	var Mtrxd3 = map[int]map[int]map[int]float64{}
	var valx3 = myfunc.Range(1, 256)
	var valy3 = myfunc.Range(1, 256)
	var valz3 = myfunc.Range(1, 64)
	for _, i := range valx3 {
		Mtrxd3[i] = make(map[int]map[int]float64)
		for _, j := range valy3 {
			Mtrxd3[i][j] = make(map[int]float64)
			for _, k := range valz3 {
				Mtrxd3[i][j][k] = 0
			}
		}
	}
	var Mtrxd4 = map[int]map[int]map[int]map[int]float64{}
	var valx4 = myfunc.Range(1, 3)
	var valy4 = myfunc.Range(1, 256)
	var valz4 = myfunc.Range(1, 256)
	var vall4 = myfunc.Range(1, 64)
	for _, i := range valx4 {
		Mtrxd4[i] = make(map[int]map[int]map[int]float64)
		for _, j := range valy4 {
			Mtrxd4[i][j] = make(map[int]map[int]float64)
			for _, k := range valz4 {
				Mtrxd4[i][j][k] = make(map[int]float64)
				for _, l := range vall4 {
					Mtrxd4[i][j][k][l] = 0
				}
			}
		}
	}

	var dummy string
	dig := []string{"0", "1", "2", "3", "4", "5", "6", "7", "8", "9"}
	partnum := StrTripArr()

	var (
		err      error
		newFile  *os.File
		newFile2 *os.File
		newFile3 *os.File
	)

	var (
		nstepx  = 256 - 1
		nstepy  = 256 - 1
		nstepz  = 64 - 1
		nstepx2 = 128 - 1
		nstepy2 = 128 - 1
		nstepz2 = 32 - 1
		pi      = math.Pi
		xcsize  = 3.125e-9
		ycsize  = 3.125e-9
		zcsize  = 3.125e-9
		ind1    = 4
		ind2    = 4
	)
	const (
		nstepxa  = 256 - 1
		nstepya  = 256 - 1
		nstepza  = 64 - 1
		nstepx2a = 128 - 1
		nstepy2a = 128 - 1
		nstepz2a = 32 - 1
	)
	var (
		gx, gy, gz, GyroSum, delta, mmgn, nrm, TurnSum, itm,
		S1, V1, S2, V2, S3, V3, S4, V4,
		FI1, FI2, FI3, FI4 float64
	)
	var (
		coorstep = Mtrxd3
		//h = Range(1, 3*(nstepxa+1)*(nstepya+1)*(nstepza+1))
		// filaa    = Mtrxd3
		// homoa    = Mtrxd3
		// turna    = Mtrxd3
		mgn  = Mtrxd4
		mgnn = Mtrxd4
		// cstep3d  = Mtrxd3
		// turn     = Mtrxd3
		fila     = Mtrxd3
		homo     = Mtrxd3
		GyroVecx = Mtrxd3
		GyroVecy = Mtrxd3
		GyroVecz = Mtrxd3
	)
	var (
	// fnamein,
	// fnameout,
	// fnameout2 string
	)
	var (
		nx0, ny0, nz0, nx, ny, nz, nnx, nny, nnz int
	)

	frag1 := `D:\mumax3.6.2\Labyr\inplane\IP\m00`
	frag5 := `D:\mumax3.6.2\Labyr\inplane\IP\gyr`
	frag4 := `D:\mumax3.6.2\Labyr\inplane\IP\fht`
	frag6 := `D:\mumax3.6.2\Labyr\inplane\IP\mr0`
	frag2 := `.ovf`
	frag3 := `.omf`

	for itm := 1; itm <= 61; itm++ {

		fnamein := frag1 + partnum[itm] + frag2
		f, err := os.Open("fnamein1")
		if err != nil {
			fmt.Println(err)
		}
		mgn = myfunc.LinePasteFile(f, err, mgn)
		f.Close()

		fnameout := frag6 + partnum[itm] + frag2
		fnameout2 := frag5 + partnum[itm] + frag2
		fnameout3 := frag4 + partnum[itm] + frag2

		os.Create(fnameout)
		newFile, err := os.OpenFile(fnameout, os.O_CREATE|os.O_APPEND|os.O_WRONLY, 0600)
		if err != nil {
			panic(err)
		}
		writeText := fmt.Sprint(`# OOMMF: rectangular mesh v1.0
# Segment count: 1
# Begin: Segment
# Begin: Header
# Title: m
# Desc:
# Desc:
# Desc:
# Desc:
# Desc:
# Desc:
# meshtype: rectangular
# meshunit: m`, "\n", `# xbase: `, xcsize, "\n",
			`# ybase: `, ycsize, "\n",
			`# zbase: `, zcsize, "\n",
			`# xstepsize: `, 2.*xcsize, "\n",
			`# ystepsize: `, 2.*ycsize, "\n",
			`# zstepsize: `, 2.*zcsize, "\n",
			`# xnodes: `, nstepx2+1, "\n",
			`# ynodes: `, nstepy2+1, "\n",
			`# znodes: `, nstepz2+1, "\n",
			`# xmin: `, 0., "\n",
			`# ymin: `, 0., "\n",
			`# zmin: `, 0., "\n",
			`# xmax: `, 2.*xcsize*(float64(nstepx2+1)), "\n",
			`# ymax: `, 2.*ycsize*(float64(nstepy2+1)), "\n",
			`# zmax: `, 2.*zcsize*(float64(nstepz2+1)), "\n",
			`# valueunit: A/m`, "\n",
			`# valuemultiplier: `, 1., "\n",
			`# ValueRangeMinMag: `, "\n",
			`# ValueRangeMaxMag: `, "\n",
			`# End: Header`, "\n",
			`# Begin: Data Text`)

		if _, err = newFile.WriteString(writeText); err != nil {
			panic(err)
		}
		for k := 1; k <= nstepz2+1; k++ {
			for j := 1; j <= nstepz2+1; j++ {
				for i := 1; i <= nstepz2+1; i++ {
					vals := fmt.Sprintf("%f %f %f", mgn[1][2*i-1][2*j-1][2*k-1], mgn[2][2*i-1][2*j-1][2*k-1], mgn[3][2*i-1][2*j-1][2*k-1])
					newFile.WriteString(vals)
				}
			}
		}
		newFile.WriteString(`# End: Data Text
		# End: Segment`)
		newFile.Close()

		for nz := 2; nz <= nstepz; nz++ {
			for ny := 2; ny <= nstepy; ny++ {
				for nx := 2; nx <= nstepx; nx++ {

					fila[nx-1][ny-1][nz-1] = mgn[3][nx-1][ny-1][nz-1]
					nx0 := 1
					ny0 := 1
					nz0 := 1
					GyroVecz[nx-1][ny-1][nz-1] = gyrox(gx, nx, ny, nz, mgn)
					GyroVecx[nx-1][ny-1][nz-1] = gyroy(gy, nx, ny, nz, mgn)
					GyroVecy[nx-1][ny-1][nz-1] = gyroz(gz, nx, ny, nz, mgn)

					GyroSum = 0
					//1 Суммируем по граням z=delta*nz и z=delta*nz0+1:
					for i := nx - 1; i <= nx; i++ {
						for j := ny - 1; j < ny; j++ {
							GyroSum -= gyroz(gz, i, j, nz-1, mgn)
							GyroSum += gyroz(gz, i, j, nz, mgn)
						}

					}

					//2 Суммируем по граням x=delta*nx0 и x=delta*nx0+1:
					for j := ny - 1; j <= ny; j++ {
						for k := nz - 1; k < nz; k++ {
							GyroSum -= gyrox(gx, nx-1, j, k, mgn)
							GyroSum += gyrox(gx, nx, j, k, mgn)
						}

					}

					//3 Суммируем по граням y=delta*ny0 и y=delta*ny:
					for i := nx - 1; i <= nx; i++ {
						for k := ny - 1; k < ny; k++ {
							GyroSum -= gyroy(gy, i, ny-1, k, mgn)
							GyroSum += gyroy(gy, i, ny, k, mgn)
						}

					}

					//Делим на 4п:
					GyroSum = GyroSum / (4. * pi)
					homo[nx-1][ny-1][nz-1] = GyroSum

				}
			}
		}

		os.Create(fnameout2)
		newFile2, err := os.OpenFile(fnameout2, os.O_CREATE|os.O_APPEND|os.O_WRONLY, 0600)
		if err != nil {
			panic(err)
		}
		writeText = fmt.Sprint(`# OOMMF: rectangular mesh v1.0
# Segment count: 1
# Begin: Segment
# Begin: Header
# Title: Gyro
# Desc:
# Desc:
# Desc:
# Desc:
# Desc:
# Desc:
# meshtype: rectangular
# meshunit: m`, "\n",
			`# xbase: `, xcsize, "\n",
			`# ybase: `, ycsize, "\n",
			`# zbase: `, zcsize, "\n",
			`# xstepsize: `, 2.*xcsize, "\n",
			`# ystepsize: `, 2.*ycsize, "\n",
			`# zstepsize: `, 2.*zcsize, "\n",
			`# xnodes: `, nstepx2, "\n",
			`# ynodes: `, nstepy2, "\n",
			`# znodes: `, nstepz2, "\n",
			`# xmin: `, 0., "\n",
			`# ymin: `, 0., "\n",
			`# zmin: `, 0., "\n",
			`# xmax: `, 2.*xcsize*(float64(nstepx2)), "\n",
			`# ymax: `, 2.*ycsize*(float64(nstepy2)), "\n",
			`# zmax: `, 2.*zcsize*(float64(nstepz2)), "\n",
			`# valueunit: A/m`, "\n",
			`# valuemultiplier: `, 1., "\n",
			`# ValueRangeMinMag: `, "\n",
			`# ValueRangeMaxMag: `, "\n",
			`# End: Header`, "\n",
			`# Begin: Data Text`)

		if _, err = newFile2.WriteString(writeText); err != nil {
			panic(err)
		}
		for k := 1; k <= nstepz2; k++ {
			for j := 1; k <= nstepz2; j++ {
				for i := 1; k <= nstepz2; i++ {
					vals := fmt.Sprintf("%f %f %f", GyroVecx[2*i-1][2*j-1][2*k-1], GyroVecy[2*i-1][2*j-1][2*k-1], GyroVecz[2*i-1][2*j-1][2*k-1])
					newFile2.WriteString(vals)
				}
			}
		}
		newFile2.WriteString(`# End: Data Text
		# End: Segment`)
		newFile2.Close()

		// 	//Закончили гомотопическое число, начинаем число вращения
		// 	for nz := 2; nz <= nstepz; nz++ { //цикл по позициям стороны [nx->nx0,nz] контура
		// 		nx0 = 1 //фиксировано: верхняя граница пленки
		// 		ny0 = 1
		// 		nz = 1
		// 		nx = nstepx + 1
		// 		ny = nstepy + 1 //меняется: положение одной из сторон прямоугольного контура

		// 		// По данным для 3D намагниченности находятся данные для искусственно вводимой 2D
		// 		// "намагниченности": проекции 3D векторов на граничную поверхность y=100нм (i=32)
		// 		// каждая нормируется на единицу. Результат записывается в массив mgnn (только для
		// 		// граничной поверхности!). Если проекция 3D вектора на границу равна нулю, 2D
		// 		// вектор также берется нулевым: сингулярность!!!).

		// 		for j := nx0; j <= nx; j++ {
		// 			for k := ny0; k < ny; k++ {
		// 				nrm = math.Sqrt(mgn[1][j][k][nz] * mgn[1][j][k][nz] * mgn[2][j][k][nz])
		// 				if nrm > 0.0000001 {
		// 					mgnn[1][j][k][nz] = mgn[1][j][k][nz] / nrm
		// 					mgnn[2][j][k][nz] = mgn[2][j][k][nz] / nrm
		// 				} else {
		// 					mgnn[1][j][k][nz] = 0
		// 					mgnn[2][j][k][nz] = 0
		// 				}
		// 			}
		// 		}
		// 		// 			Рассчитывается значение j как криволинейный интеграл по замкнутому контуру описанной
		// 		// выше формы. Роль параметра играет либо x, либо z. Произведение поля и производной
		// 		// поля аппроксимируются так:
		// 		// 	Mz(dMx/dy)(ijk)=[(1/2)(Mz(ij+1k)+Mz(ijk))]*[(Mx(ij+1k)-Mx(ijk))/(шаг сетки)].
		// 		// Далее при вычислении интегральной суммы каждое слагаемое будет содержать множитель
		// 		// (шаг сетки). Мы сразу сократим эти множители в производной и в интегральной сумме.

		// 		Turnsum := 0
		// 		for nny := 2; nny <= nstepy; nny++ {
		// 			for nnx := 2; nnx <= nstepx; nnx++ {
		// 				for j := nx - 2; j <= nx; j++ {
		// 					TurnSum += (mgnn[2][j+1][ny-1][nz]-mgnn[2][j][ny-1][nz])*
		// 						(mgnn[1][j+1][ny-1][nz]+mgnn[1][j][ny-1][nz])/2. -
		// 						(mgnn[1][j+1][ny-1][nz]-mgnn[1][j][ny-1][nz])*
		// 							(mgnn[2][j+1][ny-1][nz]+mgnn[2][j][ny-1][nz])/2.
		// 				}

		// 				for k := ny - 2; k <= ny; k++ {
		// 					TurnSum += (mgnn[2][nx][k+1][nz]-mgnn[2][nx][k][nz])*
		// 						(mgnn[1][nx][k+1][nz]+mgnn[1][nx][k][nz])/2. -
		// 						(mgnn[1][nx][k+1][nz]-mgnn[1][nx][k][nz])*
		// 							(mgnn[2][nx][k+1][nz]+mgnn[2][nx][k][nz])/2.
		// 				}

		// 				for j := nx - 2; j <= nx; j++ {
		// 					TurnSum += (mgnn[2][j+1][ny][nz]-mgnn[2][nx-1][ny][nz])*
		// 						(mgnn[1][j+1][ny][nz]+mgnn[1][nx-1][ny][nz])/2. -
		// 						(mgnn[1][j+1][ny][nz]-mgnn[1][nx-1][ny][nz])*
		// 							(mgnn[2][j+1][ny][nz]+mgnn[2][nx-1][ny][nz])/2.
		// 				}

		// 				for k := ny - 2; k <= ny; k++ {
		// 					TurnSum += (mgnn[2][nx-1][k+1][nz]-mgnn[2][nx-1][k][nz])*
		// 						(mgnn[1][nx-1][k+1][nz]+mgnn[1][nx-1][k][nz])/2. -
		// 						(mgnn[1][nx-1][k+1][nz]-mgnn[1][nx-1][k][nz])*
		// 							(mgnn[2][nx-1][k+1][nz]+mgnn[2][nx-1][k][nz])/2.
		// 				}

		// 				for nn := 2; nn <= nstepy; nn++ {
		// 					for nnx := nx - 1; nnx <= nx; nnx++ {
		// 						V1 = mgnn[1][nnx-1][nny-1][nz]*mgnn[2][nnx][nny-1][nz] -
		// 							mgnn[1][nnx][nny-1][nz]*mgnn[2][nnx-1][nny-1][nz]
		// 						S1 = mgnn[1][nnx-1][nny-1][nz]*mgnn[1][nnx][nny-1][nz] +
		// 							mgnn[2][nnx-1][nny-1][nz]*mgnn[2][nnx][nny-1][nz]
		// 						FI1 = math.Asin(S1)
		// 						TurnSum += FI1 * V1
		// 					}
		// 					V2 = mgnn[1][nnx][nny-1][nz]*mgnn[2][nnx][nny][nz] -
		// 						mgnn[2][nnx][nny-1][nz]*mgnn[1][nnx][nny][nz]
		// 					S2 = mgnn[1][nnx][nny-1][nz]*mgnn[1][nnx][nny][nz] +
		// 						mgnn[2][nnx][nny][nz]*mgnn[2][nnx][nny-1][nz]
		// 					FI2 = math.Acos(S2)
		// 					TurnSum += FI2 * V2

		// 					V3 = mgnn[1][nnx][nny][nz]*mgnn[2][nnx-1][nny][nz] -
		// 						mgnn[2][nnx][nny][nz]*mgnn[1][nnx-1][nny][nz]
		// 					S3 = mgnn[1][nnx][nny][nz]*mgnn[1][nnx-1][nny][nz] +
		// 						mgnn[2][nnx][nny][nz]*mgnn[2][nnx-1][nny][nz]
		// 					FI3 = math.Acos(S3)
		// 					TurnSum += FI3 * V3

		// 					V4 = mgnn[1][nnx-1][nny][nz]*mgnn[2][nnx-1][nny-1][nz] -
		// 						mgnn[2][nnx-1][nny][nz]*mgnn[1][nnx-1][nny-1][nz]
		// 					S4 = mgnn[1][nnx-1][nny-1][nz]*mgnn[1][nnx-1][nny][nz] +
		// 						mgnn[2][nnx-1][nny-1][nz]*mgnn[2][nnx-1][nny][nz]
		// 					FI4 = math.Acos(S4)
		// 					TurnSum += FI4 * V4
		// 					//делим на 2*pi
		// 					TurnSum /= 2. * pi
		// 					homoz[nx-1][ny-1][nz-1] = TurnSum //Значение числа вращения

		// 				}
		// 			}
		// 		}

		// 		// Разности
		// 		for nnz := 1; nnz <= nstepz-1; nnz++ {
		// 			for nny := 1; nny <= nstepy-1; nny++ {
		// 				for nnx := 1; nnx < nstepx-1; nnx++ {
		// 					fila[nnx][nny][nnz] = (filaa[nnx+1][nny][nnz] - filaa[nnx][nny][nnz]) *
		// 						(filaa[nnx][nny+1][nnz] - filaa[nnx][nny][nnz]) *
		// 						(filaa[nnx][nny][nnz+1] - filaa[nnx][nny][nnz])

		// 					homo[nnx][nny][nnz] = (homoa[nnx+1][nny][nnz] - homoa[nnx][nny][nnz]) *
		// 						(homoa[nnx][nny+1][nnz] - homoa[nnx][nny][nnz]) *
		// 						(homoa[nnx][nny][nnz+1] - homoa[nnx][nny][nnz])

		// 					turn[nnx][nny][nnz] = (turna[nnx+1][nny][nnz] - turna[nnx][nny][nnz]) *
		// 						(turna[nnx][nny+1][nnz] - turna[nnx][nny][nnz]) *
		// 						(turna[nnx][nny][nnz+1] - turna[nnx][nny][nnz])
		// 				}

		// 			}
		// 		}
		// 	}

		os.Create(fnameout3)
		newFile3, err := os.OpenFile(fnameout3, os.O_CREATE|os.O_APPEND|os.O_WRONLY, 0600)
		if err != nil {
			panic(err)
		}
		writeText = fmt.Sprint(`# OOMMF: rectangular mesh v1.0
# Segment count: 1
# Begin: Segment
# Begin: Header
# Title: FiHoTu
# Desc:
# Desc:
# Desc:
# Desc:
# Desc:
# Desc:
# meshtype: rectangular
# meshunit: m`, "\n",
			`# xbase: `, xcsize, "\n",
			`# ybase: `, ycsize, "\n",
			`# zbase: `, zcsize, "\n",
			`# xstepsize: `, 2.*xcsize, "\n",
			`# ystepsize: `, 2.*ycsize, "\n",
			`# zstepsize: `, 2.*zcsize, "\n",
			`# xnodes: `, nstepx2, "\n",
			`# ynodes: `, nstepy2, "\n",
			`# znodes: `, nstepz2, "\n",
			`# xmin: `, 0., "\n",
			`# ymin: `, 0., "\n",
			`# zmin: `, 0., "\n",
			`# xmax: `, 2.*xcsize*(float64(nstepx2)), "\n",
			`# ymax: `, 2.*ycsize*(float64(nstepy2)), "\n",
			`# zmax: `, 2.*zcsize*(float64(nstepz2)), "\n",
			`# valueunit: A/m`, "\n",
			`# valuemultiplier: `, 1., "\n",
			`# ValueRangeMinMag: `, "\n",
			`# ValueRangeMaxMag: `, "\n",
			`# End: Header`, "\n",
			`# Begin: Data Text`)

		if _, err = newFile3.WriteString(writeText); err != nil {
			panic(err)
		}
		for k := 1; k <= nstepz2; k++ {
			for j := 1; k <= nstepz2; j++ {
				for i := 1; k <= nstepz2; i++ {
					vals := fmt.Sprintf("%f %f %f", fila[2*i-1][2*j-1][2*k-1], homo[2*i-1][2*j-1][2*k-1], GyroVecz[2*i-1][2*j-1][2*k-1])
					newFile3.WriteString(vals)
				}
			}
		}
		newFile3.WriteString(`# End: Data Text
		# End: Segment`)
		newFile3.Close()

	}

	fmt.Println(math.Sqrt(2), nstepx, nstepy, nstepz, nstepx2, nstepy2, nstepz2, itm, dig, partnum[0], dummy,
		nstepxa, gx, gy, gz, GyroSum, delta, mmgn, nrm, TurnSum,
		nstepya, xcsize, ycsize, zcsize, S1, V1, S2, V2, S3, V3, S4, V4,
		nstepza, FI1, FI2, FI3, FI4,
		nstepx2a,
		nstepy2a,
		nstepz2a,
		pi, ind1, ind2, coorstep[0], frag1,
	)
}

// ErrFunc Функция ошибки
func ErrFunc(err error, text string) {
	if err != nil {
		fmt.Println(text)
	}

}

// StrTripArr массив 001 ... 999
func StrTripArr() [3000]string {
	var trip [3000]string
	for i := 0; i <= 2999; i++ {
		switch {
		case i < 10:
			trip[i] = "000" + strconv.Itoa(i)
		case i < 100 && i > 10:
			trip[i] = "00" + strconv.Itoa(i)
		case i < 1000 && i > 100:
			trip[i] = "0" + strconv.Itoa(i)
		case i > 1000:
			trip[i] = strconv.Itoa(i)
		}
	}
	return trip
}

// ФУНКЦИИ для расчета компонент плотности гиротропного вектора по
// данным о значениях намагниченности на сетке. Производные представлены
// конечно-разностными аппроксимациями; при этом квадраты шага сетки
// опущены (т.е. gx = (х-комп. гир. век.)*(шаг сетки)^2 ).
func gyrox(gx float64, i int, j int, k int, mgn map[int]map[int]map[int]map[int]float64) float64 {
	if math.Abs(mgn[3][i][j][k]) > 0.57 {
		gx = ((mgn[2][i][j+1][k] - mgn[2][i][j][k]) * (mgn[1][i][j][k+1] - mgn[1][i][j][k])) -
			(mgn[2][i][j][k+1]-mgn[2][i][j][k])*(mgn[1][i][j+1][k]-mgn[1][i][j][k])/
				mgn[3][i][j][k]

	} else if math.Abs(mgn[1][i][j][k]) > 0.57 {
		gx = ((mgn[3][i][j+1][k] - mgn[3][i][j][k]) * (mgn[2][i][j][k+1] - mgn[2][i][j][k])) -
			(mgn[3][i][j][k+1]-mgn[3][i][j][k])*(mgn[2][i][j+1][k]-mgn[2][i][j][k])/
				mgn[1][i][j][k]
	} else {
		gx = ((mgn[1][i][j+1][k] - mgn[1][i][j][k]) * (mgn[3][i][j][k+1] - mgn[3][i][j][k])) -
			(mgn[1][i][j][k+1]-mgn[1][i][j][k])*(mgn[3][i][j+1][k]-mgn[3][i][j][k])/
				mgn[2][i][j][k]
	}
	return gx
}

func gyroy(gy float64, i int, j int, k int, mgn map[int]map[int]map[int]map[int]float64) float64 {
	if math.Abs(mgn[3][i][j][k]) > 0.57 {
		gy = ((mgn[2][i][j][k+1] - mgn[2][i][j][k]) * (mgn[1][i+1][j][k] - mgn[1][i][j][k])) -
			(mgn[2][i+1][j][k]-mgn[2][i][j][k])*(mgn[1][i][j][k+1]-mgn[1][i][j][k])/
				mgn[3][i][j][k]

	} else if math.Abs(mgn[1][i][j][k]) > 0.57 {
		gy = ((mgn[3][i][j][k+1] - mgn[3][i][j][k]) * (mgn[2][i+1][j][k] - mgn[2][i][j][k])) -
			(mgn[3][i+1][j][k]-mgn[3][i][j][k])*(mgn[2][i][j][k+1]-mgn[2][i][j][k])/
				mgn[1][i][j][k]
	} else {
		gy = ((mgn[1][i][j][k+1] - mgn[1][i][j][k]) * (mgn[3][i+1][j][k] - mgn[3][i][j][k])) -
			(mgn[1][i+1][j][k]-mgn[1][i][j][k])*(mgn[3][i][j][k+1]-mgn[3][i][j][k])/
				mgn[2][i][j][k]
	}
	return gy
}

func gyroz(gz float64, i int, j int, k int, mgn map[int]map[int]map[int]map[int]float64) float64 {
	if math.Abs(mgn[3][i][j][k]) > 0.57 {
		gz = ((mgn[2][i+1][j][k] - mgn[2][i][j][k]) * (mgn[1][i][j+1][k] - mgn[1][i][j][k])) -
			(mgn[2][i][j+1][k]-mgn[2][i][j][k])*(mgn[1][i+1][j][k]-mgn[1][i][j][k])/
				mgn[3][i][j][k]

	} else if math.Abs(mgn[1][i][j][k]) > 0.57 {
		gz = ((mgn[3][i][j][k+1] - mgn[3][i][j][k]) * (mgn[2][i+1][j][k] - mgn[2][i][j][k])) -
			(mgn[3][i+1][j][k]-mgn[3][i][j][k])*(mgn[2][i][j][k+1]-mgn[2][i][j][k])/
				mgn[1][i][j][k]
	} else {
		gz = ((mgn[1][i+1][j][k] - mgn[1][i][j][k]) * (mgn[3][i][j+1][k] - mgn[3][i][j][k])) -
			(mgn[1][i][j+1][k+1]-mgn[1][i][j][k])*(mgn[3][i+1][j+1][k]-mgn[3][i][j][k])/
				mgn[2][i][j][k]
	}
	return gz
}