Theoretical Physics II: Soft Matter - Heinrich-Heine-Universität Düsseldorf

Institut für Theoretische Physik II: Weiche Materie - HHU Düsseldorf

Fortgeschrittenen-Praktikum: Weiche Materie - Archiv - Milch und Mayo Lösung (Programm)

Quelltext des Programms:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<iostream.h>
#include<time.h>

#define Z 50 //Anzahl der Abstandsintervalle//

double a, //Halbachse//
       b,//Halbachse//
       beta, //Parameter der Oberflaechenspannung//
       F, // Gesamtflaeche einer Ellipse //
       b_1,b_2,b_1alt,b_2alt, //Halbachsen auf der Verbindungslinie der Schwerpunkte//
       d,d_alt, // Abstand der Schwerpunkte der Ellipsen//
       d_max, //max Abstand der Schwerpunkte der Ellipsen//
       r_k[Z+1], //Gitter fuer Abstand der Schwerpunkte der Ellipsen//
       U_1,U_2, //Umfaenge der beiden Ellipsen//
       E_mitt[Z]={0}, //mittlere Energie//
       E, //neue Energie//
       E_alt, //Energie//
       E_ges[Z]={0}, //Gesamtenergie im Abstand d//
       var_dmax, //Prozentsatz um den d maximal variiert wird//
       var_bmax; //Prozentsatz um den b maximal variiert wird//

long double pw,pt; //Variablen fuer Zufallszahlen//

long int t; //Zaehler in main//

long int N, //Anzahl der Variationen im Abstand d //
         h_d[Z]={0}; //Moeglichkeiten im Abstand d//

double Umfang(double);
void Eingabe(void);
double zufall(void);
void start(void);

double zufall(void) //generiert Zufallszahl zwischen 0 und 1//
{
long double pm;

pm=rand();
pm=pm/RAND_MAX;
return pm;
}

void Eingabe(void)
{
char Taste;
int i;

   cout << "Moechten Sie die Flaeche (F) oder die grossen Halbachsen der";
   cout << " Ellipsen (h) vorgeben? " ;
   cout << " Druecken Sie bitte die entsprechenden Tasten: ";
   cin>>Taste;
   if ((Taste=='f') || (Taste=='F'))
    {
    cout<<" Geben Sie die Flaeche an: ";
    cin>>F;
    }
   else
    {
     cout<<" Geben Sie die Halbachse a an: ";
     cin>>a;
     cout<<" Geben Sie die Halbachse b an: ";
     cin>>b;
     F=M_PI*a*b;
    }
   cout<<" Geben Sie den Parameter Beta an: ";
   cin>>beta;
   cout<<" Geben Sie die Anzahl der Variationen der Parameter b_1,b_2 und r an: ";
   cin>>N;
   cout<<" Geben Sie den maximalen Abstand der Schwerpunkte ein: ";
   cin>>d_max;
   cout<<" Um wieviel Prozent des maximalen Schwerpunktsabstands sollen die";
   cout<<" Halbachsen der Ellipsen maximal variieren?";
   cin>>var_bmax;
   cout<<" Bitte warten";
   for (i=0;i<=Z;i++)
    {
     r_k[i]=i*d_max/Z;
    }
   }

void start(void) //Bestimmt die Ausgangskonfiguration zufaellig//
{
do
{
d_alt=zufall()*d_max;
}
while (d_alt==0);
b_1alt=b_2alt=sqrt(F/M_PI);
if (b_1alt+b_2alt>d_alt) b_1alt=b_2alt=d_alt/2;
E_alt=beta*(Umfang(b_1alt)+Umfang(b_2alt));
}

double var_gr(double alt,double var_max,double d_max) //Variation der Groessen//
{
double gr;

gr=alt+(zufall()-0.5)*(2*var_max/100)*d_max;

if (gr<=0) gr=alt+zufall()*(var_max/100)*d_max;
return gr;
}

void ueberlap(double c_1,double c_2,double b_1alt,double b_2alt,double ab,double d_alt)
{
if ((c_1+c_2)>ab)
{
   b_1=b_1alt;
   b_2=b_2alt;
   d=d_alt;
}
}

double Umfang(double b) //Bestimmung des Umfangs der Ellipse//
{
double lambda, //Ellipsenparameter//
U; //Umfang der Ellipsen//

a=F/(M_PI*b); //Bestimmung der zweiten Halbachse//
lambda=(fabs(a-b))/(a+b);
U=M_PI*(a+b)*(64-3*pow(lambda,4))/(64-16*pow(lambda,2));
return U;
}

double E_warsch(double En,double U_1,double U_2)
//Monte Carlo Schritte//
{
long double p, //Hilsvariable//

E=beta*(U_1+U_2);
p=exp(En-E);
return p;
}

double Entscheid(long double p,long double pt,double gr,double alt)
{
double gr_hilf;

if (p>=pt) gr_hilf=gr;
else gr_hilf=alt;
return gr_hilf;
}

void Ener_prue(double d,double Ener) //Bestimmt die Gesamtenergie im Abstand d//
          //und Bestimmt die Haufigkeit im Abstand d//
{
int i;
for (i=0;i<=Z-1;i++)
{
   if ((d>r_k[i]) && (d<=r_k[i+1]))
   {
    E_ges[i]=E_ges[i]+Ener;
    h_d[i]=h_d[i]+1;
   }
}
}

void Ener_mit(void)
{
int i;

for (i=0;i<=Z-1;i++)
{
   if (h_d[i]==0)
   {
    E_mitt[i]=0;
   }
   else E_mitt[i]=E_ges[i]/h_d[i];
}
}

void speichern(void)
{
FILE *datei_1,*datei_2,*datei_3,*datei_4;
double s;

int i;
datei_1=fopen("Ellipse.dat","wt");
datei_2=fopen("Energie.dat","wt");
datei_3=fopen("lnEllip.dat","wt");
datei_4=fopen("Umfang.dat","wt");
for (i=0;i<=Z-1;i++)
   {
    s=(r_k[i]+r_k[i+1])/2;
    fprintf(datei_1,"%le %li ",s,h_d[i]);
    if (E_mitt[i]!=0)
    {
     fprintf(datei_2,"%le %le ",s,E_mitt[i]);
     fprintf(datei_4,"%le %le ",s,E_mitt[i]/beta);
    }
    if ((h_d[i])!=0)
    {
     fprintf(datei_3,"%le %le ",s,-log(h_d[i]));
    }
   }
fclose(datei_1);
fclose(datei_2);
fclose(datei_3);
fclose(datei_4);
}

double abstand(void)
{
double a;
do
{
   a=zufall()*d_max;
}
while (a==0);
return a;
}

void main(void)
{
Eingabe();
start();
Ener_prue(d,E_alt);
for (t=1;t<=N;t++)
{
   d=abstand();
   b_1=var_gr(b_1alt,var_bmax,d_max);
   b_2=var_gr(b_2alt,var_bmax,d_max);
   ueberlap(b_1,b_2,b_1alt,b_2alt,d,d_alt);
   U_1=Umfang(b_1);
   U_2=Umfang(b_2);
   pw=E_warsch(E_alt,U_1,U_2);
   pt=zufall();
   b_1alt=Entscheid(pw,pt,b_1,b_1alt);
   b_2alt=Entscheid(pw,pt,b_2,b_2alt);
   d_alt=Entscheid(pw,pt,d,d_alt);
   E_alt=beta*(Umfang(b_1alt)+Umfang(b_2alt));
   Ener_prue(d_alt,E_alt);
}
Ener_mit();
speichern();
cout<<" Fertig ";
}

zurück