Added gradient calculation

isotables/isotopapp/views.py

@@ -1,14 +1,18 @@
+import os
 from math import pi
 from io import BytesIO
 import base64
+
+from bokeh.models.axes import LinearAxis
 from django.shortcuts import render
 from django.utils.safestring import mark_safe
 import re
 
+
+import numpy as np
 from bokeh.plotting import figure
 from bokeh.embed import components
-from bokeh.models import Label, Node, MathML
-
+from bokeh.models import Label, Node, MathML, Range1d, Span
 
 from isotopapp.models import Isotope
 # Create your views here.
@@ -17,6 +21,10 @@ def home(request):
     isotopes = [i for i in Isotope.objects.all() if (i.gamma != 0 or i.stable)]
     return render(request, 'home.html', {'isotopes': [[f"{i.n_nucleons}{i.symbol}", mark_safe(f"<sup>{i.n_nucleons}</sup>{i.symbol}")] for i in isotopes],})
 
+def sfg(request):
+    isotopes = [i for i in Isotope.objects.all() if (i.gamma != 0 or i.stable)]
+    return render(request, 'sfg.html', {'isotopes': [[f"{i.n_nucleons}{i.symbol}", mark_safe(f"<sup>{i.n_nucleons}</sup>{i.symbol}")] for i in isotopes],})
+
 def extract_isotope_parts(isotope_str):
     """Extracts the number and element from an isotope string (e.g., '23Na')."""
     print(isotope_str)
@@ -131,3 +139,96 @@ def result(request):
     else:
         ans = []
     return render(request, 'result.html', {'ans': ans, 'script': script, 'div': div})
+
+def position(request):
+    print(request)
+    n1, element1 = extract_isotope_parts(request.GET.get('isotope1'))
+    isotope1 = Isotope.objects.filter(symbol=element1, n_nucleons=n1).get()
+    print(os.path.abspath("."))
+    data = np.loadtxt(request.GET.get('magnet'))
+    _z_coords = data[:,0]
+    _fields = data[:,1]
+    _gradients = data[:,2]
+
+    freq = float(request.GET.get('freq'))
+    field_T = freq / isotope1.gamma
+
+    sample_diameter = 5e-3
+    gradient = float(request.GET.get('gradient'))
+    # create a plot (bokeh)
+
+    TOOLTIPS = [
+        ("index", "$index"),
+        ("(x,y)", "($x, $y)"),
+        ("desc", "@desc"),
+    ]
+    plot = figure(outer_width=500, outer_height=500, match_aspect=False, tooltips=TOOLTIPS)
+    frame_left = Node(target="frame", symbol="left", offset=5)
+    frame_bottom = Node(target="frame", symbol="bottom", offset=0)
+    frame_top = Node(target="frame", symbol="top", offset=0)
+    #plot.ellipse(x=[0], y=[0], width=5, height=5, color="#D5D9FF", alpha=0.8, line_width=1, line_color="black")
+    #plot.ellipse(x=[0], y= [2.5], width=5, height=5, color="#D5D9FF", alpha=0.4, line_width=1, line_color="black")
+    #plot.ellipse(x=[0], y=[-2.5], width=5, height=5, color="#D5D9FF", alpha=0.4, line_width=1, line_color="black")
+    plot.xaxis[0].axis_label = 'z in mm'
+    plot.yaxis[0].axis_label = 'B0 in T'
+    plot.line(x=_z_coords*1e3, y=_fields, color="navy", line_width=2)
+    plot.yaxis[0].axis_label_text_color = "navy"
+    plot.y_range = Range1d(-10, 10)
+
+
+    plot.extra_y_ranges['gradient'] = Range1d(-200, 200)
+    plot.line(x=_z_coords * 1e3, y=_gradients, color="crimson", line_width=2 , y_range_name="gradient")
+    ax2 = LinearAxis(y_range_name="gradient", axis_label="g in T/m")
+    ax2.axis_label_text_color = "crimson"
+    plot.add_layout(ax2, 'left')
+
+    # find fields where _fields == field:
+    pos_threshold_indices = np.where(np.diff(np.sign(_fields + field_T)) != 0)[0]
+    neg_threshold_indices = np.where(np.diff(np.sign(_fields - field_T)) != 0)[0]
+    indices = np.concatenate((pos_threshold_indices, neg_threshold_indices))
+
+
+    for p in indices:
+        vline = Span(location=_z_coords[p]*1e3, dimension='height', line_color='gray', line_width=1)
+        plot.renderers.extend([vline])
+        label = Label(x=_z_coords[p]*1e3, y=_fields[p]+1, text=f"{_z_coords[p]*1e3:.2f},{_gradients[p]:.1f}T/m", text_baseline="middle",
+                text_align="center", text_font_size="11pt")
+        plot.add_layout(label)
+    plot.circle(x=_z_coords[indices]*1e3, y=_fields[indices], radius=2, color="gray")
+    plot.circle(x=_z_coords[indices]*1e3, y=_gradients[indices], radius=2, color="gray",y_range_name="gradient" )
+    hline1 = Span(location=field_T, dimension='width', line_color='gray', line_width=1)
+    hline2 = Span(location=-field_T, dimension='width', line_color='gray', line_width=1)
+    plot.renderers.extend([hline1, hline2])
+
+    close_isotopes = []
+    for isotope in Isotope.objects.all():
+        if isotope.gamma == 0: continue
+        if not isotope.stable: continue
+        z = (freq/isotope.gamma-field_T)/gradient
+        if abs(z) <= sample_diameter:
+            i_info = isotope_info(isotope, field_T)
+            #i_info[3] = f"{z*1e3:.1f} mm"
+            close_isotopes.append(i_info)
+            #plot.rect(x=[0], y=[z*1e3], width=5, height=0.2, color="black", alpha=0.6)
+            #label = Label(x=2.6, y=z*1e3, text=f"{isotope.n_nucleons}{isotope.symbol}", text_baseline="middle", text_align="left", text_font_size="16pt")
+            #plot.add_layout(label)
+
+
+    citation = Label(
+        x=frame_left,
+        y=frame_bottom,
+        anchor="bottom_left",
+        text=f"{isotope1.n_nucleons}{isotope1.symbol}: {freq:.1f} MHz\ng={gradient:.1f} T/m\n5 mm sample dia.",
+        padding=5,
+        border_radius=5,
+        border_line_color="#D5D9FF", background_fill_color="white",
+    )
+    plot.add_layout(citation)
+
+    # boke plot
+    script, div = components(plot)
+    ans = sorted(close_isotopes, key=lambda x: float(x[3]))
+
+    return render(request, 'posiiton.html', {'ans': ans, 'script': script, 'div': div})
+
+