diff --git a/test.cpp b/test.cpp
index dbc1c24..1c89d36 100644
--- a/test.cpp
+++ b/test.cpp
@@ -37,8 +37,12 @@ static vector<double> energies;
 static vector<string> orbString;
 
 
+static vector<string> orbitalStr;
+static vector<double> AList;
+static vector<vector<double>> enLevelList;
+
 GLFWwindow* g_window;
-ImVec4 clear_color = ImVec4(0.45f, 0.55f, 0.60f, 1.00f);
+ImVec4 clear_color = ImVec4(0.427, 0.702f, 0.243f, 1.00f);
 int g_width;
 int g_height;
 
@@ -92,48 +96,92 @@ void loop(){
     ImGui::End();
   }
 
-  // // // Demo Window
-  // if (true){
-  //   ImGui::SetNextWindowPos(ImVec2(650, 20), ImGuiCond_FirstUseEver); // Normally user code doesn't need/want to call this because positions are saved in .ini file anyway. Here we just want to make the demo initial state a bit more friendly!
-  //   ImGui::ShowDemoWindow();
-  // }
-
   //  Woods-Saxon window
   {
-    ImGui::SetNextWindowSize(ImVec2(1000, 1000), ImGuiCond_FirstUseEver);
-    static float V0 = -45, R0 = 3.5, a0 = 0.6;
-    static float VSO = 28, RSO = 3.5, aSO = 0.6;
+    ImGui::SetNextWindowSize(ImVec2(600, 1000), ImGuiCond_FirstUseEver);
+    static float V0 = -45, r0 = 1.25, R0 = 3.5, a0 = 0.6;
+    static float VSO = 28, rSO = 1.25, RSO = 3.5, aSO = 0.6;
     static int Z = 0, nStep = 300;
-    static float Rc = 3.5, dr = 0.1;
-
+    static float Rc = 3.5, rc = 1.25, dr = 0.1;
+    static int A = 20;
+    static float mass = 939.565;
 
     ImGui::Begin("Woods-Saxon Calculation");
 
-    ImGui::SliderFloat("V0 [MeV]", &V0, -100, 0);
-    ImGui::SliderFloat("R0 [fm]", &R0, 1, 10);
-    ImGui::SliderFloat("a0 [fm]", &a0, 0.1, 2);
-    ImGui::SliderFloat("VSO [MeV]", &VSO, 0, 40);
-    ImGui::SliderFloat("RSO [fm]", &RSO, 1, 10);
-    ImGui::SliderFloat("aSO [fm]", &aSO, 0.1, 2);
-    ImGui::SliderInt("Z ", &Z, 0, 100);
-    ImGui::SliderFloat("Rc [fm] ", &Rc, 1, 10);
-    ImGui::SliderInt("nStep", &nStep, 100, 400);
-    ImGui::SliderFloat("dr [fm]", &dr, 0.01, 1);
+    static int e = 0;
+    ImGui::RadioButton("Arbitary", &e, 0); ImGui::SameLine();
+    ImGui::RadioButton("Isotope", &e, 1);
+    ImGui::Separator();
 
+    ImGui::PushItemWidth(100);
+    static int nucleon = 0;
+    ImGui::RadioButton("Neutron", &nucleon, 0); ImGui::SameLine();
+    ImGui::RadioButton("Proton", &nucleon, 1);
+    
+    ImGui::BeginDisabled(e==0); ImGui::DragInt("A ", &A, 1, 1, 200); ImGui::EndDisabled();
+    ImGui::SameLine();
+    ImGui::DragInt("Z", &Z, 1, 0, 100);
+
+    ImGui::BeginDisabled(e==1); ImGui::DragFloat("Rc [fm] ", &Rc, 0.01, 1, 10); ImGui::EndDisabled();
+    ImGui::SameLine();
+    ImGui::BeginDisabled(e==0); 
+    ImGui::DragFloat("rc [fm] ", &rc, 0.01, 1, 10); 
+    ImGui::EndDisabled();
+
+    if( e == 0 ) {A = 1; rc = Rc / pow(A, 1./3);}
+    if( e == 1 ) {
+      rc = 1.25;
+      r0 = 1.25;
+      rSO = 1.25;
+      Rc = rc * pow(A, 1./3);
+    }
+
+    ImGui::Text("Eff. mass : %.3f ", mass);
+
+    ImGui::Separator();
+
+    ImGui::DragFloat("V0 [MeV]", &V0, 0.1, -100, 0); 
+
+    ImGui::BeginDisabled(e==1); ImGui::DragFloat("R0 [fm] ", &R0, 0.01, 1, 10); ImGui::EndDisabled();
+    ImGui::SameLine();
+    ImGui::BeginDisabled(e==0); ImGui::DragFloat("r0 [fm] ", &r0, 0.01, 1, 10); ImGui::EndDisabled();
+    if( e == 0 ) {A = 1; r0 = R0 / pow(A, 1./3);}
+    if( e == 1 ) R0 = r0 * pow(A, 1./3);
+
+    ImGui::DragFloat("a0 [fm]", &a0, 0.01, 0.1, 2); 
+    ImGui::DragFloat("VSO [MeV]", &VSO, 0.1, 0, 40);
+    ImGui::BeginDisabled(e==1); ImGui::DragFloat("RSO [fm] ", &RSO, 0.01, 1, 10); ImGui::EndDisabled();
+    ImGui::SameLine();
+    ImGui::BeginDisabled(e==0); ImGui::DragFloat("rSO [fm] ", &rSO, 0.01, 1, 10); ImGui::EndDisabled();
+    if( e == 0 ) {A = 1; rSO = RSO / pow(A, 1./3);}
+    if( e == 1 ) RSO = rSO * pow(A, 1./3);
+    ImGui::DragFloat("aSO [fm]", &aSO, 0.01, 0.1, 2);
+
+    ImGui::Separator();
+    ImGui::DragInt("nStep", &nStep, 50, 100, 400); 
+    ImGui::SameLine();
+    ImGui::DragFloat("dr [fm]", &dr, 0.001, 0.01, 0.1);
 
     ImGui::PushStyleVar(ImGuiStyleVar_FramePadding, ImVec2(30, 10));
-    if( ImGui::Button("Cal WS Level") ){
+    if( ImGui::Button("Cal WS Levels") ){
       WoodsSaxon ws;
 
-      if( Z == 0 ) {
-        ws.SetNucleus(1,1);
+      if( nucleon == 0) {
         ws.IsNeutron();
       }else{
-        ws.SetNucleus(1, Z);
         ws.IsProton();
-        ws.SetRc(Rc);
       }
 
+      if( e == 0) {
+        ws.SetNucleus(1,Z);
+      }else{
+        ws.SetNucleus(A,Z);
+        mass = 931.5 * (ws.GetA() )/(1.008664 + ws.GetA());
+        ws.SetMass(mass);
+      }
+
+      if( Z > 0 ) ws.SetRc(Rc);
+
       ws.SetV0( V0 );
       ws.SetR0( R0 );
       ws.Seta0( a0 );
@@ -173,22 +221,22 @@ void loop(){
     }
     ImGui::PopStyleVar();
 
+    ImGui::SetNextItemOpen(true, ImGuiCond_Always);
     if( ImGui::TreeNode("WS Energy Levels:") ){
       for( int i = 0; i < energies.size(); i++){
-
         char buf[32];
         sprintf(buf, "%24.12f MeV %s", energies[i],orbString[i].c_str());
         if( ImGui::Selectable(buf, selected == i) ) selected = i;
-        //ImGui::Text("%12.5f MeV %s", ws.energy[i], ws.orbString[i].c_str());
       }
       ImGui::TreePop();
     }
 
-    if( ImPlot::BeginPlot("Plot") ){
+    if( ImPlot::BeginPlot("Plot", ImVec2(ImGui::GetWindowWidth()-20, 400)) ){
       ImPlot::SetupLegend(ImPlotLocation_SouthEast);
       ImPlot::SetupAxes("r [fm]"," Energy [MeV]");
-      ImPlot::SetupAxesLimits(0, 10, floor(V0*1.1), 1);
+      ImPlot::SetupAxesLimits(0, 10, floor(V0*1.1), 10);
       ImPlot::SetupAxisLimitsConstraints(ImAxis_X1, 0, 30);
+      ImPlot::SetupAxisLimitsConstraints(ImAxis_Y1, floor(V0*1.1), abs(floor(V0*1.1)));
       ImPlot::PlotLine("Central", xValues, WSCValues, nPt);
       ImPlot::PlotLine("S-O", xValues, WSSOValues, nPt);
 
@@ -207,6 +255,235 @@ void loop(){
     ImGui::End();
   }
 
+  //*============= WS Trend
+  {
+    ImGui::SetNextWindowSize(ImVec2(600, 800), ImGuiCond_FirstUseEver);
+    ImGui::SetNextWindowPos(ImVec2(100, 200), ImGuiCond_FirstUseEver);
+    static float V0 = -45, r0 = 1.25, a0 = 0.6;
+    static float VSO = 28, rSO = 1.25,  aSO = 0.6;
+    static int Z = 0, N = 1, nStep = 300;
+    static float rc = 1.25, dr = 0.1;
+    static int A[2] = {10, 20};
+    static int NRange[2] = {1, 20};
+    static int ZRange[2] = {1, 20};
+    static float mass = 939.565;
+    static float kappa = -1;
+
+    ImGui::Begin("Woods-Saxon Calculation (Range)");
+
+    static int e = 0;
+    ImGui::RadioButton("Stablility", &e, 0); ImGui::SameLine();
+    ImGui::RadioButton("var. A", &e, 1); ImGui::SameLine();
+    ImGui::RadioButton("var. N", &e, 2); ImGui::SameLine();
+    ImGui::RadioButton("var. Z", &e, 3);
+    ImGui::Separator();
+
+    ImGui::PushItemWidth(150);
+    static int nucleon = 0; //neutron or proton
+    if( e != 1 ){
+      ImGui::RadioButton("Neutron", &nucleon, 0); ImGui::SameLine();
+      ImGui::RadioButton("Proton", &nucleon, 1);
+    }else{
+      ImGui::Text("Neutron only");
+    }
+
+    if( e == 0 || e == 1 ){
+      ImGui::DragInt2("A range", A, 1, 1, 300);
+    }
+
+    if( e == 2){
+
+      ImGui::DragInt("Z", &Z, 0, 0, 100);
+      ImGui::DragInt2("N range", NRange, 1, 1, 300);
+    }
+
+    if( e == 3){
+      ImGui::DragInt("N", &N, 1, 1, 200);
+      ImGui::DragInt2("Z range", ZRange, 1, 0, 300);
+    }
+
+    if( e != 1){
+      ImGui::DragFloat("kappa", &kappa, 0.1, -10, 10); 
+    }
+
+    ImGui::Separator();
+
+    ImGui::DragFloat("V0 [MeV]", &V0, 0.1, -100, 0); 
+    ImGui::DragFloat("r0 [fm] ", &r0, 0.01, 1, 10); 
+    ImGui::DragFloat("a0 [fm]", &a0, 0.01, 0.1, 2); 
+    ImGui::DragFloat("VSO [MeV]", &VSO, 0.1, 0, 40);
+    ImGui::DragFloat("rSO [fm] ", &rSO, 0.01, 1, 10); 
+    ImGui::DragFloat("aSO [fm]", &aSO, 0.01, 0.1, 2);
+
+    ImGui::Separator();
+    ImGui::DragInt("nStep", &nStep, 50, 100, 400); 
+    ImGui::SameLine();
+    ImGui::DragFloat("dr [fm]", &dr, 0.001, 0.01, 0.1);
+
+
+    int ARange[2] ;
+    int fixedNZ;
+    if( e == 0 || e == 1 ) {ARange[0] = A[0]; ARange[1] = A[1]; fixedNZ = 0;}
+    if( e == 2 ) {ARange[0] = NRange[0] + Z; ARange[1] = NRange[1] + Z; fixedNZ = Z;}
+    if( e == 3 ) {ARange[0] = ZRange[0] + N; ARange[1] = ZRange[1] + N; fixedNZ = N;}
+    int maxA = ARange[1];
+    if( maxA < A[0]) {
+      maxA = ARange[0];
+      ARange[0] = ARange[1];
+      ARange[1] = maxA;
+    }
+
+    orbitalStr.clear();
+    orbitalStr.push_back("0s1/2");   ///2
+    orbitalStr.push_back("0p3/2");
+    orbitalStr.push_back("0p1/2");   ///8
+    if( maxA > 8 ) {
+      orbitalStr.push_back("0d5/2");
+      orbitalStr.push_back("1s1/2");
+      orbitalStr.push_back("0d3/2");   ///20
+    }
+    if( maxA > 20 ){
+      orbitalStr.push_back("0f7/2");   ///28
+      orbitalStr.push_back("1p3/2");
+      orbitalStr.push_back("0f5/2");
+    }
+    if( maxA > 40 ){
+      orbitalStr.push_back("1p1/2");   ///40
+      orbitalStr.push_back("0g9/2");   ///50
+      orbitalStr.push_back("0g7/2");   
+      orbitalStr.push_back("1d5/2");   
+      orbitalStr.push_back("0h11/2"); 
+      orbitalStr.push_back("1d3/2");   
+    }
+    if( maxA > 82 ){
+      orbitalStr.push_back("2s1/2");   ///82
+      orbitalStr.push_back("0h9/2"); 
+      orbitalStr.push_back("1f7/2"); 
+      orbitalStr.push_back("0i13/2"); 
+      orbitalStr.push_back("2p3/2"); 
+      orbitalStr.push_back("1f5/2"); 
+    }
+    if( maxA > 126 ){
+      orbitalStr.push_back("2p1/2");    ///126 
+      orbitalStr.push_back("1g9/2"); 
+      orbitalStr.push_back("0i11/2"); 
+      orbitalStr.push_back("0j15/2"); 
+      orbitalStr.push_back("0j15/2"); 
+      orbitalStr.push_back("2d5/2"); 
+      orbitalStr.push_back("3s1/2"); 
+      orbitalStr.push_back("2d3/2"); 
+      orbitalStr.push_back("1g7/2");    ///184 
+    }
+
+    static float progress = 0;
+    static bool cal = false;
+    static bool finsihed = false;
+    static int ANZ = ARange[0];
+
+    if( ImGui::Button("Cal WS Trend") ){
+      cal = true;
+      finsihed = false;
+      AList.clear();
+      enLevelList.clear();
+      ANZ = ARange[0];
+    }
+
+    if(cal) {
+      WoodsSaxon ws;
+
+      ws.Setr0(r0);
+      ws.SetrSO(rSO);
+      ws.Setrc(rc);
+      ws.SetVSO(VSO);
+      ws.SetRange2(0.001, dr, nStep);
+      if( nucleon == 0 ) {
+        ws.IsNeutron();
+      }else{
+        ws.IsProton();
+      }
+      
+      if( ANZ <= ARange[1] ){
+         
+        AList.push_back(ANZ);
+        if( e == 0){
+          ws.SetNucleus ( ANZ, 0.008 + 0.495 * ANZ - 0.00076 * ANZ *ANZ + 1.4e-6 * ANZ * ANZ * ANZ); // aproximate stable vallage
+          ws.CalRadius();
+          ws.SetV0 ( V0 * ( 1 +  kappa * ( 2.0* ws.GetN() - ws.GetA() ) / ws.GetA() ) ); 
+        }
+
+        if( e == 1){
+          ws.SetNucleus ( ANZ, 0);
+          ws.CalRadius();
+          ws.SetV0( V0 );
+        }
+
+        if( e == 2 ){ /// change in N, ANZ = N
+          ws.SetNucleus( ANZ , fixedNZ );
+          ws.CalRadius();
+
+          ws.SetV0 ( V0 * ( 1 +  kappa * ( 2.0* ws.GetN() - ws.GetA() ) / ws.GetA() ) ); 
+        }
+    
+        if( e == 3 ){ /// change in Z, ANZ = Z
+          ws.SetNucleus( ANZ , ANZ - fixedNZ );
+          ws.CalRadius();
+          ws.SetV0 ( V0 * ( 1 +  kappa * ( 2.0* ws.GetN() - ws.GetA() ) / ws.GetA() ) ); 
+        }
+
+        double reducedMass = 931.5 * (ws.GetA() )/(1.008664 + ws.GetA());
+        ws.SetMass(reducedMass);
+
+        ws.ClearVector();
+        ws.CalWSEnergies(false, 7, 500, 1e-7, 400, 0.2, false);
+
+        //ws.PrintWSParas();
+        // ws.PrintEnergyLevels();
+
+        vector<double> exList;
+        for( int i = 0; i < orbitalStr.size(); i++){
+          bool isMatched = false;
+          for( int j = 0; j < (ws.energy).size(); j++){
+            if( orbitalStr[i] == ws.orbString[j] ) {
+              exList.push_back(ws.energy[j]);
+              isMatched = true;
+            }
+          }
+          if( isMatched == false ) exList.push_back(0);
+        }
+        enLevelList.push_back(exList);
+
+        progress = (ANZ - ARange[0]) * 1.0 / (ARange[1] - ARange[0]);
+
+        ANZ ++;
+      }else{
+        cal = false;
+        finsihed = true;
+      }
+    }
+    ImGui::SameLine(); ImGui::ProgressBar(progress);
+    
+    if( finsihed ){
+      if( ImPlot::BeginPlot("Plot", ImVec2(ImGui::GetWindowWidth()-20, 400)) ){
+        ImPlot::SetupLegend(ImPlotLocation_East, ImPlotLegendFlags_Outside);
+        ImPlot::SetupAxes("A"," Energy [MeV]");
+        ImPlot::SetupAxesLimits(ARange[0], ARange[1], -50, 0);
+        // ImPlot::SetupAxisLimitsConstraints(ImAxis_X1, ARange[0], ARange[1]);
+        // ImPlot::SetupAxisLimitsConstraints(ImAxis_Y1, floor(V0)*2, 0);
+        const double * kaka = AList.data();
+        for( int i = 0; i < orbitalStr.size(); i++){
+          vector<double> jaja;
+          for( int j = 0; j < AList.size(); j++) jaja.push_back(enLevelList[j][i]);
+          const double * haha = jaja.data();
+          ImPlot::PlotLine(orbitalStr[i].c_str(), kaka, haha, AList.size());
+        }
+
+        ImPlot::EndPlot();
+      }
+    }
+
+    ImGui::End();
+  }
+
   ImGui::Render();
 
   int display_w, display_h;