# Objective: linear
# Constraints: convex quadratic
# Feasible set: convex

################################################################
# Narrow band 3-dimensional beam pattern optimization
################################################################

param pi := 4*atan(1);

param d := 0.5;
set Coords dimen 3; data coords.dat;
set Els := setof {(i,j,k) in Coords} (i*d, j*d, k);
set Eli := setof {(i,j,k) in Coords} (k);
display card(Els);
display card(Eli);

param width := sqrt(card(Els));
#param delta := 1/(2*width); display delta;
param delta := 1/(width); display delta;
param delmax := ceil(2/(sqrt(3)*delta));
set P := setof {i in -delmax..delmax, j in -delmax..delmax:
		(i + j/2)^2 + 3*j^2/4 <= 1/delta^2} 
		(delta*(i+j/2), delta*j*sqrt(3)/2);

# The x,y,z components of the look direction
param lx := sin(pi*22.5/180)*cos(pi*30/180);  #22.5 degrees that's half way
param ly := sin(pi*22.5/180)*sin(pi*30/180);  #22.5 degrees
#param lx := sin(pi*0/180)*cos(pi*30/180);  #0 degrees
#param ly := sin(pi*0/180)*sin(pi*30/180);  #0 degrees
param lz := sqrt(1-lx^2-ly^2);

param beamwidth := 7;  # in degrees (actually halfwidth i.e. from center)
param cosbeamwidth := cos(beamwidth*pi/180);

# (px,py,pz) is inward pointing, whereas 
# (lx,ly,lz) is outward pointing
# p dot -l should be <= a threshold
set P_Sidelobes := setof {(px,py) in P: 
		-px*lx - py*ly + sqrt(abs(1-px^2-py^2))*lz <= cosbeamwidth}
		(px,py); 

printf {(px, py) in P}: "%6.2f %6.2f \n", px, py > p.out;
printf {(elx, ely, eli) in Els}: "%6.2f  %6.2f %d \n", elx, ely, eli > el.out;
display card(P);

var w_real {Eli};
var w_imag {Eli};

var A_real {(px,py) in P union {(-lx,-ly)} } = 
    sum {(elx,ely,eli) in Els} (w_real[eli]*cos(-2*pi*(px*elx + py*ely))
			      - w_imag[eli]*sin(-2*pi*(px*elx + py*ely))
			       );

var A_imag {(px,py) in P union {(-lx,-ly)} } = 
    sum {(elx,ely,eli) in Els} (w_real[eli]*sin(-2*pi*(px*elx + py*ely))
			      + w_imag[eli]*cos(-2*pi*(px*elx + py*ely))
			       );

#sine  [px,py,elx,ely,eli]       

param rho := 10^(-15/10);   # 15 dB

minimize L_two:
    sum {(px,py) in P_Sidelobes} 
        ((A_real[px,py]^2 + A_imag[px,py]^2)/sqrt(1.01-(px^2+py^2)));

#subject to chebychev {(px,py) in P_Sidelobes}:
#    A_real[px,py]^2 + A_imag[px,py]^2 <= rho;

subject to gain_real:
    A_real[-lx,-ly] = 1;

subject to gain_imag:
    A_imag[-lx,-ly] = 0;

let {eli in Eli} w_real[eli] := Uniform01();
solve;

printf {eli in Eli}: "%d %6.2f %6.2f \n",
        eli, w_real[eli], w_imag[eli] > nbsup.out;

# this finer grid is just for finer detail in the output
param delta2 := delta/3;
param delmax2 := ceil(2/(sqrt(3)*delta2));
set P2 := setof {i in -delmax2..delmax2, j in -delmax2..delmax2:
		(i + j/2)^2 + 3*j^2/4 <= 1/delta2^2} 
		(delta2*(i+j/2), delta2*j*sqrt(3)/2);

param A2_real {P2};
let {(px,py) in P2} A2_real[px,py] := 
    sum {(elx,ely,eli) in Els} (w_real[eli]*cos(-2*pi*(px*elx + py*ely))
			      - w_imag[eli]*sin(-2*pi*(px*elx + py*ely))
			       );

param A2_imag {P2};
let {(px,py) in P2} A2_imag[px,py] := 
    sum {(elx,ely,eli) in Els} (w_real[eli]*sin(-2*pi*(px*elx + py*ely))
			      + w_imag[eli]*cos(-2*pi*(px*elx + py*ely))
			       );
#param r;
#param x {0..1,0..1};
#param y {0..1,0..1};
#for {i in -delmax2-1..delmax2+1, j in -delmax2-1..delmax2+1} {
#    for {ii in 0..1, jj in 0..1} {
#        let x[ii,jj] := delta2*((i+ii)+(j+jj)/2);
#        let y[ii,jj] := delta2*(j+jj)*sqrt(3)/2;
#    }
#    if (x[0,0]^2 + y[0,0]^2 <= 1 &&
#        x[0,1]^2 + y[0,1]^2 <= 1 &&
#        x[1,0]^2 + y[1,0]^2 <= 1) then {
#
#        for {ii in 0..1, jj in 0..1: ii+jj<=1} {
#            let r := max(10*log10(A2_real[x[ii,jj],y[ii,jj]]^2
#                                + A2_imag[x[ii,jj],y[ii,jj]]^2) + 50, 0);
#            printf "%6.2f %6.2f %6.2f ", x[ii,jj], y[ii,jj], r > faces;
#        }
#        printf "\n" > faces;
#    };
#    if (x[1,1]^2 + y[1,1]^2 <= 1 &&
#        x[0,1]^2 + y[0,1]^2 <= 1 &&
#        x[1,0]^2 + y[1,0]^2 <= 1) then {
#
#        for {ii in 0..1, jj in 0..1: ii+jj>=1} {
#            let r := max(10*log10(A2_real[x[ii,jj],y[ii,jj]]^2
#                                + A2_imag[x[ii,jj],y[ii,jj]]^2) + 50, 0);
#            printf "%6.2f %6.2f %6.2f ", x[ii,jj], y[ii,jj], r > faces;
#        }
#        printf "\n" > faces;
#    };
#}

set IJ := setof {i in -delmax2..delmax2, j in -delmax2..delmax2:
		(i + j/2)^2 + 3*j^2/4 <= 1/delta2^2 ||
		(i+1 + j/2)^2 + 3*j^2/4 <= 1/delta2^2 ||
		(i-1 + j/2)^2 + 3*j^2/4 <= 1/delta2^2 ||
		(i + (j+1)/2)^2 + 3*(j+1)^2/4 <= 1/delta2^2 ||
		(i + (j-1)/2)^2 + 3*(j-1)^2/4 <= 1/delta2^2 
		} 
		(i,j);

set IJ0 := setof {i in -delmax2..delmax2, j in -delmax2..delmax2:
		(i + j/2)^2 + 3*j^2/4 <= 1/delta2^2 
		} 
		(i,j);

param x {(i,j) in IJ} := delta2*(i+j/2);
param y {(i,j) in IJ} := delta2*j*sqrt(3)/2;

printf {(i,j) in IJ0: x[i+1,j]^2+y[i+1,j]^2<=1 && x[i,j+1]^2+y[i,j+1]^2<=1}:
 "%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n",
 x[i,j],   y[i,j],   max(10*log10(A2_real[x[i,  j],  y[i,  j  ]]^2+
		                  A2_imag[x[i,  j],  y[i,  j  ]]^2)+50,0),
 x[i+1,j], y[i+1,j], max(10*log10(A2_real[x[i+1,j],  y[i+1,j  ]]^2+
		                  A2_imag[x[i+1,j],  y[i+1,j  ]]^2)+50,0),
 x[i,j+1], y[i,j+1], max(10*log10(A2_real[x[i,  j+1],y[i,  j+1]]^2+
		                  A2_imag[x[i,  j+1],y[i,  j+1]]^2)+50,0)
 > faces;

printf {(i,j) in IJ0: x[i-1,j]^2+y[i-1,j]^2<=1 && x[i,j-1]^2+y[i,j-1]^2<=1}:
 "%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n",
 x[i,j],   y[i,j],   max(10*log10(A2_real[x[i,  j],  y[i,  j  ]]^2+
		                  A2_imag[x[i,  j],  y[i,  j  ]]^2)+50,0),
 x[i-1,j], y[i-1,j], max(10*log10(A2_real[x[i-1,j],  y[i-1,j  ]]^2+
		                  A2_imag[x[i-1,j],  y[i-1,j  ]]^2)+50,0),
 x[i,j-1], y[i,j-1], max(10*log10(A2_real[x[i,  j-1],y[i,  j-1]]^2+
		                  A2_imag[x[i,  j-1],y[i,  j-1]]^2)+50,0)
 > faces;