add scattering files

riemann.py

@@ -1,22 +1,96 @@
-# This script demonstrates the prime numbers are quasicrystalline.
+# This script demonstrates the distribution of prime numbers is quasicrystalline.
+
+# A quasicrystal is an arrangement of atoms without strict periodicity, but having a point-like diffraction spectrum
+
+# An infinite periodic lattice   
 
 # prime counting function, pi(n), returns the number of primes less than or equal to n.
 # so, pi(0) = 0, pi(1) = 1, pi(6) = 4,
 
+from numpy import log as ln
+from mpmath import li
+from scipy.fft import fft
+import math
+
 def is_prime(n):
-    # could also use int(sqrt(n))+1 as the upper limit:w
+    if n < 2: return False
 
-    for i in range(n, int(n/2)+1):
-        if i % n == 0:
-            return 0
-    return 1
+    for i in range(2,int(math.sqrt(n))+1):
+        if (n%i) == 0:
+            return False
+    return True
 
 def quick_pi(n):
-    return sum([int(is_prime(i)) for i in range(n)]) 
-
-if __name__ == "__main__":
-    print(quick_pi(1))
-    print(quick_pi(1000))
-    print(quick_pi(100000))
-    print(quick_pi(1003))
-
+    return sum([int(is_prime(i)) for i in range(1,n)]) 
+
+def is_primes(n):
+    return [int(is_prime(i)) for i in range(0,n)]
+
+def primes(n):
+    # returns all primes up to n
+    return [i for i in range(1,n) if is_prime(i)]
+
+def ln_transform(seq):
+    return [s/ln(s) for s in seq]
+
+def li_transform(seq):
+    return [s/li(s) for s in seq]
+
+def li_transform2(seq):
+    return [i/li(i) if s !=0 else 0 for i, s in enumerate(seq)]
+
+def RZF_trunc(x, N):
+    return sum([1/(s^x) for s in range(N)])
+
+def RZF_zeros(N, tol=1E-6):
+    # return the first N RZF zeros
+    zeros = []
+    while len(zeros) < N:
+        # get next zero
+        pass         
+    return zeros
+
+
+#for z in [100, 1000, 10000, 100000, 1000000]:
+#    p = primes(z)
+
+#print(primes(100))    
+#print(ln_transform(primes(100)))
+#print(len(ln_transform(primes(100)))/ln_transform(primes(100))[-1])
+#print(len(ln_transform(primes(1000)))/ln_transform(primes(1000))[-1])
+#print(len(ln_transform(primes(10000)))/ln_transform(primes(10000))[-1])
+#print(len(ln_transform(primes(100000)))/ln_transform(primes(100000))[-1])
+#print(len(ln_transform(primes(1000000)))/ln_transform(primes(1000000))[-1])
+#print(len(ln_transform(primes(10000000)))/ln_transform(primes(10000000))[-1])
+
+
+#print(li(10))
+#input()
+#print('done')
+from scipy.special import zeta
+#print(zeta(10))
+#input()
+
+# Consider the integral from (0,inf) of the 
+X = 500    
+import matplotlib.pyplot as plt
+#plt.plot([i for i in range(X)],is_primes(X), label="primes")
+print(li_transform2(is_primes(X)))
+
+
+
+
+#plt.plot([li_transform(primes(X))], [int(bool(primes(X))])) 
+
+
+#plt.plot([i/ln(i) for i in range(X*1)],primes(X*1), label="log transformed primes")
+#plt.plot([i for i in range(X)], [quick_pi(i)/X for i in range(X)], label="pi(x)")
+#plt.plot([i/ln(i) for i in range(X*1)], [quick_pi(i)/(X*1) for i in range(X*1)], label="log transformed pi(x)")
+#plt.plot(fft(primes(3000)))
+
+#plt.legend()
+
+plt.show()
+
+
+

scattering.f

@@ -0,0 +1,76 @@
+program scart
+        implicit real*8(a-h,o-z)
+        dimension nprim(1000),freal(1000),fimag(1000)
+        dimension kpoint(100),npsht(1000,500),nprimh(5000)
+
+        open(8,file='prime',status='unknown',form='formatted')
+        open(9,file='cosida',status='unknown',form='formatted')
+c
+c  this program feeds kpoint in unit of 2*pi/a'
+c  calculate cos(k*nprim) and sine(k*nprim)
+c
+        read(8,100)(nprimh(i),i=1,61)
+  100   format(5x,10i10/(5x,10i10))
+        do 77 ipr = 1,61
+        nprim(ipr) = nprimh(ipr)
+  77    continue
+        write(9,100)(nprim(i),i=1,61)
+c
+c  shift the prime numbers to the average
+c  the cos ans sin are calculated roughly symmetric
+c
+        pi = 3.14159
+        twopi = 2*pi
+
+        idsum = 0
+        limsi = 61
+        lacli = 53
+        do 1 isum = 1,limsi
+        idsum = idsum + nprim(isum)
+  1     continue
+        nshift = idsum/limsi
+        do 6 imovc = 1,3
+        if(imovc.eq.0) then
+        nmovd = 0
+        else
+        nmovd = nprimh(imovc)
+        endif
+        mshitol = nshift + nmovd
+        do 4 ish = 1,lacli
+        npsht(ish,imovc) = nprim(ish) - nshift
+  4     continue
+        do 5 ish = 1,lacli
+        nprim(ish) = npsht(ish,imovc)
+  5     continue
+        write(9,104)idsum,imove,nshift,nmovd,mshitol
+ 104    format(1x,'6 idsum,imoven,shift,nmovd,mshitol',2x,5i10)
+        write(9,102)(npsht(ish,imovc),ish=1,50)
+ 102    format(5x,10i10/(5x,10i10))
+c
+c  do the scattering
+c
+
+       do 2  ik = 1,20
+       gk = twopi*ik*0.1
+       sumco = 0.d0
+       sumim = 0.0d0
+       summag = 0.0d0       do 3 ipo = 1,lacli
+       pprim = nprim(ipo)
+       phase = gk*pprim
+       sumco = sumco + dcos(phase)
+       sumsi = sumsi + dsin(phase)
+       if(ipo.le.5) then
+       write(9,103)nprim(ipo),phase,sumco,sumsi
+  103  format(1x,'nprim,phase,sumco,sumsi',2x,i10,3f15.5)
+       endif
+       sumcosq = sumco**2
+       sumsisq = sumsi**2
+       sumcosi = sumcosq + sumsisq
+       summag = summag + dsqrt(sumcosi)
+   3   continue
+       write(9,101)ik,gk,sumco,sumsi,summag
+  101  format(1x,'k,gk,sumco,sumsi,summag',2x,i10,4f12.5)
+   2   continue
+   6   continue
+       stop
+       end

scattering.py

@@ -0,0 +1,81 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Define arrays
+prime_numbers = np.zeros(1000)  # Prime numbers
+shifted_prime_numbers = np.zeros((1000,500))  # Shifted prime numbers
+prime_number_holder = np.zeros(5000)  # Holder for prime numbers read from file
+
+# Open the input files
+with open('prime', 'r') as f_prime, open('cosida', 'w') as f_cosida:
+    # Read prime numbers from 'prime' file
+    prime_number_holder[:61] = np.fromfile(f_prime, sep=" ", count=61)
+
+    # Assign prime numbers
+    prime_numbers[:61] = prime_number_holder[:61]
+
+    # Write prime numbers to 'cosida' file
+    np.savetxt(f_cosida, prime_numbers[:61], fmt='%10i')
+
+    # Constants
+    pi = 3.14159
+    two_pi = 2*pi
+
+    sum_id = 0
+    limit_prime = 61
+    limit_array = 53
+    sum_id = np.sum(prime_numbers[:limit_prime])
+    shift_number = sum_id/limit_prime
+
+    # Scattering amplitude data for plot
+    scattering_amplitude = []
+    momentum = []
+
+    for move_count in range(1, 4):
+        move_number = prime_number_holder[move_count] if move_count != 0 else 0
+        shift_total = shift_number + move_number
+
+        # Shift the prime numbers
+        shifted_prime_numbers[:limit_array,move_count] = prime_numbers[:limit_array] - shift_number
+        prime_numbers[:limit_array] = shifted_prime_numbers[:limit_array,move_count]
+
+        f_cosida.write(f'{sum_id:<10}{move_count:<10}{shift_number:<10}{move_number:<10}{shift_total:<10}\n')
+        np.savetxt(f_cosida, shifted_prime_numbers[:50, move_count], fmt='%10i')
+
+        # Scattering process
+        for wave_number in range(1, 21):
+            wave_vector = two_pi * wave_number * 0.1  # wave vector (k) in units of 2*pi/a
+            sum_cosine = 0  # Sum of cosines
+            sum_sine = 0  # Sum of sines
+            sum_magnitude = 0  # Sum of magnitudes
+
+            for prime_index in range(1, limit_array+1):
+                prime = prime_numbers[prime_index]
+                phase = wave_vector * prime  # Phase = k * prime number
+                sum_cosine += np.cos(phase)  # Summing cosines
+                sum_sine += np.sin(phase)  # Summing sines
+
+                if prime_index <= 5:
+                    f_cosida.write(f'{prime_numbers[prime_index]:<10}{phase:<15.5f}{sum_cosine:<15.5f}{sum_sine:<15.5f}\n')
+
+                sum_cosine_squared = sum_cosine ** 2
+                sum_sine_squared = sum_sine ** 2
+                sum_cosine_sine = sum_cosine_squared + sum_sine_squared
+                sum_magnitude += np.sqrt(sum_cosine_sine)  # Magnitude of the scattering amplitude
+
+            f_cosida.write(f'{wave_number:<10}{wave_vector:<12.5f}{sum_cosine:<12.5f}{sum_sine:<12.5f}{sum_magnitude:<12.5f}\n')
+
+            # Append computed data for plotting
+            scattering_amplitude.append(sum_magnitude)
+            momentum.append(wave_vector)
+
+# After the computations, we can now plot the scattering amplitude
+plt.figure(figsize=(10, 6))
+plt.plot(momentum, scattering_amplitude, label='Scattering Amplitude')
+plt.xlabel('Momentum (k in units of 2*pi/a)')
+plt.ylabel('Scattering Amplitude')
+plt.title('Scattering Amplitude vs Momentum')
+plt.grid(True)
+plt.legend()
+plt.show()
+