while 1:
    try:
        shell_number = int(input("Enter shell number for estimation (100): "))
    except:
        print("Shell number set at 100.")
        shell_number = 100

    init = [0, 0, 1]
    triplets = [init]

    while (triplets[-1][0 and 1 and 2])<shell_number:
        lasttrip=(triplets[-1])
        if lasttrip[0]==lasttrip[1]:
            if lasttrip[1]==lasttrip[2]:
                triplets.append([0, 0, lasttrip[2]+1])
            elif lasttrip[1]<lasttrip[2]:
                triplets.append([0, lasttrip[1]+1, lasttrip[2]])
        elif lasttrip[0]<lasttrip[1] and lasttrip[1]<=lasttrip[2]:
            triplets.append([lasttrip[0]+1, lasttrip[1], lasttrip[2]])

    #Generates an equal length list of distances from the reference atom (0,0,0) to each triplet.      

    distset=[]
    for i in range(len(triplets)):
        coordset=triplets[i]
        distance = ((coordset[0]**2)+(coordset[1]**2)+(coordset[2]**2))**(1/2)
        distset.append(distance)

    #print("List of triplets used:", triplets)

    #Separates into odd (Cl) which adds and even (Na) which subtracts

    denomlist = []
    for i in range(len(triplets)):
        coordset = triplets[i]
        if coordset[0] == coordset[1]:
            if coordset[1] == coordset[2]:
                denom = 8
            else:
                denom = 2
        elif coordset[1] == coordset[2]:
            denom = 4
        else:
            denom = 2
        denomlist.append(denom)

    #print("List of denominators used for atom based on fraction in shell (for faster convergence):", tripletdenoms)

    eqdistset = []
    for i in range(len(triplets)):
        coordset=triplets[i]
        if coordset[0] == 0:
            if coordset[1] == 0:
                eqdist = 6
            elif coordset[1] == coordset[2]:
                eqdist = 12
            else:
                eqdist = 24
        elif coordset[0] == coordset[1]:
            if coordset[1] == coordset[2]:
                eqdist = 8
            else:
                eqdist = 24
        else:
            if coordset[1] == coordset[2]:
                eqdist = 24
            else:
                eqdist = 48
        eqdistset.append(eqdist)

    #print("List of N atoms at the same distance:", eqdistset)

    addlist = []
    for i in range(len(triplets)):
        coordset=triplets[i]
        if (coordset[0]+coordset[1]+coordset[2])%2 == 1:
            if max(coordset) == shell_number:
                addlist.append(eqdistset[i]/(denomlist[i]*distset[i]))
            else:
                addlist.append(eqdistset[i]/distset[i])
        else:
            if max(coordset) == shell_number:
                addlist.append(-eqdistset[i]/(denomlist[i]*distset[i]))
            else:
                addlist.append(-eqdistset[i]/distset[i])


    #print("List of final constants used for adding: %.4f" % elem in addlist for i in range(len(addlist)))
    print("Final calculated Madelung constant for NaCl is:", sum(addlist))