subroutine script_fermi_lines(scriptunit, seedname, fermi_energy_list)
!================================================!
use w90_constants, only: dp
implicit none
real(kind=dp), allocatable, intent(in) :: fermi_energy_list(:)
integer, intent(in) :: scriptunit
character(len=50), intent(in) :: seedname
write (scriptunit, '(a)') &
"# Energy level for isocontours (typically the Fermi level)"
write (scriptunit, '(a,f12.6)') "ef=", fermi_energy_list(1)
write (scriptunit, '(a)') " "
write (scriptunit, '(a)') &
"bands=np.loadtxt('"//trim(seedname)//"-kslice-bands.dat')"
write (scriptunit, '(a)') "numbands=bands.size//num_pt"
write (scriptunit, '(a)') "if square:"
write (scriptunit, '(a)') &
" bbands=bands.reshape((dimy,dimx,numbands))"
write (scriptunit, '(a)') " for i in range(numbands):"
write (scriptunit, '(a)') " Z=bbands[:,:,i]"
write (scriptunit, '(a)') " pl.contour(x_coord,y_coord,Z," &
//"[ef],colors='black')"
write (scriptunit, '(a)') "else:"
write (scriptunit, '(a)') " bbands=bands.reshape((num_pt," &
//"numbands))"
write (scriptunit, '(a)') " bandint=[]"
write (scriptunit, '(a)') " grid_x, grid_y = np.meshgrid(xint,yint)"
write (scriptunit, '(a)') " for i in range(numbands):"
write (scriptunit, '(a)') " bandint.append(interpolate.griddata" &
//"((points_x,points_y), bbands[:,i], (grid_x,grid_y), " &
//"method='nearest'))"
write (scriptunit, '(a)') " pl.contour(grid_x,grid_y," &
//"bandint[i],[ef],colors='black')"
end subroutine script_fermi_lines