' temhorn_lens

Sub Main ()
    Dim B1, B2, H1, H2, L1, L2 As Double
    B1 = RestoreDoubleParameter("B1")
    B2 = RestoreDoubleParameter("B2")
    H1 = RestoreDoubleParameter("H1")
    H2 = RestoreDoubleParameter("H2")
    L1 = RestoreDoubleParameter("L1")
    L2 = RestoreDoubleParameter("L2")

    Dim nx, ny, nz As Integer
    nx = RestoreDoubleParameter("nx")
    ny = RestoreDoubleParameter("ny")
    nz = RestoreDoubleParameter("nz")

    Dim epsMax As Double
    epsMax = RestoreDoubleParameter("epsmax")

    Dim var As Integer
	var = RestoreDoubleParameter("lensvariant")

    Dim lh, lb, l As Double
	lh = H1 * L2 / (H2 - H1)
	lb = B1 * L2 / (B2 - B1)
	If lb < lh Then
		l = lb
	Else
		l = lh
	End If


	Dim x0, y0, z0 As Double
	x0 = -B2/2
	y0 = -H2/2

	If var = 0 Then
		z0 = -lh
	Else
		z0 = -l
	End If


	Dim dx, dy, dz As Double
	dx = B2 / nx
	dy = H2 / ny
	dz = (-z0 + L2) / nz

	Dim ix, iy, iz As Integer

	For ix = 0 To nx-1
	For iy = 0 To ny-1
	For iz = 0 To nz-1
		Dim x, y, z As Double
		x = x0 + (ix+0.5)*dx
		y = y0 + (iy+0.5)*dy
		z = z0 + (iz+0.5)*dz

		Dim minX, maxX, minY, maxY, minZ As Double
		minX = x - dx/2
		maxX = x + dx/2
		minY = y - dy/2
		maxY = y + dy/2
		minZ = z - dz/2

		If (2*( maxY - H1/2) / (H2 - H1) < minZ / L2) And _
		   (2*( maxX - B1/2) / (B2 - B1) < minZ / L2) And _
		   (2*(-minY - H1/2) / (H2 - H1) < minZ / L2) And _
		   (2*(-minX - B1/2) / (B2 - B1) < minZ / L2) Then

			Dim eps As Double

			If var = 0 Then
				eps = epsMax * (z+l)^2 / (x*x + y*y + (z+l)^2)
			ElseIf var = 1 Then
				Dim R1sq, R2sq As Double
				R1sq = (L2 + lb)^2 + (B2/2)^2
				R2sq = (L2 + lh)^2 + (H2/2)^2

				Dim cosAlphaSq, cosBetaSq As Double
				cosAlphaSq = (lb+L2)^2 / R1sq
				cosBetaSq = (lh+L2)^2 / R2sq

				Dim tgPhiSq, cosXiSq As Double
				tgPhiSq = y^2 / (lh+z)^2
				cosXiSq = (lb+z)^2 / ((lb+z)^2 + x^2)

				eps = epsMax * R1sq*cosAlphaSq / (R1sq*cosAlphaSq/cosXiSq + R2sq*cosBetaSq*tgPhiSq)
			ElseIf var = 2 Then
				Dim cosPhiSq, tgXiSq As Double
				cosPhiSq = (lh+L2)^2 / R2sq
				tgXiSq = x^2 /(lb+z)^2

				eps = epsMax * R2sq*cosBetaSq / (R1sq*cosAlphaSq*tgXiSq + R2sq*cosBetaSq/cosPhiSq)

			End If



			eps = Fix(eps*50)/50


Dim epsName As String
epsName = "" & eps
epsName = Replace(epsName, ",", ".")


			' Create material
			With Material
	        	.Reset
	        	.Folder   "temhorn_lens"
	        	.Name     "material_eps=" & epsName
	        	.FrqType  "all"
	        	.Type     "Normal"
	        	.Epsilon  eps
	        	.Mue      1
	        	.Create
	    	End With

			' Create brick
			With Brick
	            .Reset
	            .Component "temhorn_lens"
	            .Name      "solid_" & ix & "_" & iy & "_" & iz
	            .Material  "temhorn_lens/material_eps=" & epsName
	            .Xrange x0 + ix*dx, x0 + (ix+1)*dx
	            .Yrange y0 + iy*dy, y0 + (iy+1)*dy
	            .Zrange z0 + iz*dz, z0 + (iz+1)*dz
	            .Create
	        End With
        End If
	Next
	Next
	Next

End Sub