Unit Cells:
Rhombohedral ↔ Hexagonal

On[Assert]

Transformations

P[a]=2/3a+1/3b+1/3c
P[b]=-1/3a+1/3b+1/3c
P[c]=-1/3a-2/3b+1/3c

P[a]==a1&&P[b]==b1&&P[c]==c1;
sol=Solve[%,{a,b,c}][[1]];
{Q[a1],Q[b1],Q[c1]}={a,b,c}/.sol

From hexagonal to rhombohedral

RH={
a->1/2{1,Sqrt[3],0},
b->1/2{1,-Sqrt[3],0},
c->{0,0,h}
}
Assert[a.b==-1/2&&a.c==0&&b.c==0/.RH]

RR={
a1->P[a],
b1->P[b],
c1->P[c]
}/.RH
Rall=Join[RH,RR];
Assert[a1.a1==b1.b1==c1.c1/.Rall]

Rhombehedral unit cell angle as a function of the height of the hexagonal cell

a1.b1/a1.a1/b1.b1/.Rall//Simplify

OO={0,0,0};
Graphics3D[{
{AbsoluteThickness[3],
{
Arrow[{{OO,a},{OO,b}}],
{Darker[Blue],Arrow[{OO,c}]},
Line[{{a,2a+b},{b,2b+a},{2a+b,2a+2b},{2b+a,2b+2a}}]
},
{Red,
Arrow[{{OO,a1},{OO,b1},{OO,c1}}],Line[{{a1,a1+b1},{a1,a1+c1},{b1,b1+a1},{b1,b1+c1},{c1,c1+a1},{c1,c1+b1}}],
Line[{a1+b1,a1+b1+c1,b1+c1}],
Line[{a1+c1,a1+b1+c1}]
},
{Red,Dashed,
Line[{a1,a1-b1}],
Line[{b1,b1-c1}],
Line[{a1,a1-c1,2a1-c1,2a1-2c1}],
Line[{a1-c1,a1+b1-c1,a1+b1-2c1}],
Line[{a1+b1-c1,b1-c1}],
Line[{2a1-c1,2a1-b1-c1}],
Line[{a1,2a1,2a1-c1}],
Line[{2a1,2a1-b1,2a1-b1-c1}],
Line[{2a1-b1,a1-b1}],
Line[{2a1-c1,2a1-c1+b1,2a1-2c1+b1,a1-2c1+b1}],
Line[{2a1-2c1+b1,2a1-2c1}],
Line[{2a1-c1+b1,a1-c1+b1,a1+b1}]
}
}
},
Boxed->False,PlotLabel->Style["Hexagonal cell in black and blue\nRhombohedral cell in plain and red\nLines of same colour have the same length",18]]/.Rall/.h->3/2

Assert[Norm[a-a1]==Norm[a1]/.Rall]

From rhombohedral to hexagonal

RR1={
a1->{Cos[θ] Sqrt[3]/2,-Cos[θ]/2,Sin[θ]},
b1->{0,Cos[θ],Sin[θ]},
c1->{-Cos[θ] Sqrt[3]/2,-Cos[θ]/2,Sin[θ]}
}
Assert[Simplify[a1.a1==1&&b1.b1==1&&c1.c1==1&&a1.b1==b1.c1==c1.a1/.RR1,Assumptions:>θ ϵ Reals]]

Rhombohedral unit cell angle as a function of the elevation angle shared by the 3 sides

a1.b1/a1.a1/b1.b1/.RR1//Simplify

RH1={
a->Q[a1],
b->Q[b1],
c->Q[c1]
}/.RR1
Assert[a.a==b.b/.RH1]
Assert[a.c==0/.RH1]
Assert[b.c==0/.RH1]
Assert[Simplify[a.b/Norm[a]/Norm[b]==-1/2/.RH1,Assumptions:>0<θ<Pi/2]]