23 Mart 2012 Cuma

Çanak Anten Formülü


 
Bir önceki yazıda ispatladığımız gibi yer-sabit (geostationary) uydular ekvator düzleminde yer alan 42.164 km yarıçaplı bir çember üzerinde dolanırlar. Şu anda bu yörüngede 300'den fazla uydu vardır. Bu uyduların konumları ekvator üzerine izdüşümlerinin meridyeni ile verilir. Örnek olarak yer-sabit bir uydu olan TÜRKSAT uydusu 42 derece doğu meridyeni üzerinde yer alır. Peki dünya üzerinde herhangi bir noktadan herhangi bir uyduya çevirdiğimiz çanak anteninin yönü ne olmalıdır? Vektör işlemleri ve koordinat sistemi dönüşümleri ile bu soru da kolayca cevaplandırabilir.

Dünyanın mükemmel bir küre olduğunu varsayalım; (bu varsayımın hassasiyetine en son bakarız.)



İki koordinat sistemi seçiyoruz. Birincisi orjini dünyanın merkezinde olan, z-ekseni kuzey kutbundan çıkan, x-ekseni 0 derece meridyeni ile ekvatorun kesişim noktasından çıkan kartezyen koordinat sistemi olsun. Bu koordinat sisteminde yer-sabit uydunun konumunu vektörü ile; dünya üzerinde paraleli ve meridyeninde bulunan gözlemcinin konumunu ise vektörü ile gösterelim. Bu durumda uyduya bakan gözlemci vektörü yönünde bakmaktadır.

Seçtiğimiz ikinci koordinat sistemi gözlemcinin konum vektörü için küresel koordinat sistemidir. Bu sistem aynı zamanda bu gözlemcinin yerel koordinat sistemidir. Öyle ki vektörü gök kubbenin "doruğunu" (zenith) gösterirken vektörü coğrafi kuzeyi (manyetik kuzeyle karışmasın) vektörü ise coğrafi doğu istikametini göstermektedir. Şapka sembolü birim vektörleri ifade etmektedir.

Vektörlerimizi bu koordinat sistemlerinin birim vektörleri cinsinden yazalım.






Burada R bir önceki yazıda hesapladığımız 42.164 km. olup r ise dünyanın yarıçapı olan 6371 km. dir. (Küre kabulunde bu ortalama rakamı alıyoruz.) Bunun haricinde tanımladığımız küresel birim vektörlerinin kartezyen birim vektörler cinsinden ifadelerine de ihtiyacımız olacak:








Esas yapmamız gereken vektörünün küresel sistemdeki birim vektörler yönünde izdüşümlerini bulmakdır. Böylece hangi coğrafi yönde ve ufuktan kaç derece yukarı bakılması gerektiğini bulmuş oluruz.

Vektörler arasındaki açıyı bulmak için iki vektörü skaler çarpım yaptırıp boylarına böleceğiz. Bu işlem bize aranılan açının kosinüsünü verir.

Önce doruk ile yapılan açıya bakalım. Bu açıya ismini verelim. Bu durumda






eşitliğinin sağındaki ifadeyi hesaplamamız gerekmektedir. Bu ifadenin ters kosinüsünü alıp 90 dereceden çıkarırsak gözlemcinin uyduyu ufuktan kaç derece yukarıda (elevation) araması gerektiğini buluruz.

Peki gözlemcinin yüzünü hangi yöne çevirmesi gerekmektedir?  Bunu bulmak için vektörünün (kuzey) ile (doğu) üzerinde izdüşümlerinin oranını hesaplayacağız. Bu oran bize yönü bildiren açısının tanjantını verecektir.



(Burada beta eğer 0 derece ise kuzeyi, 90 ise doğuyu, 180 ise güneyi, 270 ise batıyı ifade eder.)

Vektörlerin hepsinin kartezyen koordinatların birim vektörleri cinsinden ifadelerini yerine koymak kalıyor. Eğer sözüme güvenirseniz sonuçta çıkan formüllerin hiç de sevimli ve akılda kalıcı şeyler olmadığını söyleyebilirim. Dolayısı ile daha pratik olan şeyi tercih ediyorum ve bu hesabı bizim için yapan bir bilgisayar kodu veriyorum. Her zamanki gibi python dilinde. Sonuçlar sadece bu işi yapan bir web sitesi olan www.dishpointer.com ile tamamen uyumludur. Dolayısı ile küre varsayımımızın iyi işlediği söylenebilir. 

NOT: Aynı programda R vektörü dünya üzerinde bir noktayı gösterecek şekilde değiştirilir, omega ve alfa iptal edilirse program bu noktaya doğru bakan yönü bulmak için de pekala kullanılabilir. (ör: Kıble yönü vs...)

---------------
import math as m  # Trigonometri icin lazim

def carp(a,b):           # Once vektorleri skaler carpmayi ogretelim bilgisayara
    return a[0]*b[0]+a[1]*b[1]+a[2]*b[2]

def skaler(c,a):         # Simdi de bir vektoru bir skaler ile carpmayi
    d = [ ]
    d.append(c*a[0])
    d.append(c*a[1])
    d.append(c*a[2])
    return d

def topla(a,b):         # Iki vektoru toplamayi da ogretelim.
    c = [ ]
    c.append(a[0]+b[0])
    c.append(a[1]+b[1])
    c.append(a[2]+b[2])
    return c

def boy(a):              # Son olarak bir vektorun boyu nasil hesaplanir gosterelim
    return (a[0]*a[0]+a[1]*a[1]+a[2]*a[2])**0.5

r_boy = 6371.0             # Dunya yaricapi
R_boy = 42164.0          # Uydu mesafesi

print 'KUZEY PARALELLERI POZITIF GUNEY PARALELLERI NEGATIF GIRILMELIDIR'

teta = input('Paralel: ')*m.pi/180    # Girilen aciyi dereceden radyana cevirir

print 'DOGU MERIDYENLERI POZITIF BATI MERIDYENLERI NEGATIF GIRILMELIDIR'

fi = input('Meridyen: ')*m.pi/180

print 'DOGU MERIDYENLERI POZITIF BATI MERIDYENLERI NEGATIF GIRILMELIDIR'

omega = input('Uydu meridyeni: ')*m.pi/180

# Birim vektorleri tanitalim

r_bir = [m.cos(teta)*m.cos(fi),m.cos(teta)*m.sin(fi),m.sin(teta)]
teta_bir = [-1*m.sin(teta)*m.cos(fi),-1*m.sin(teta)*m.sin(fi),m.cos(teta)]
fi_bir = [-1*m.sin(fi),m.cos(fi),0]


# Uydunun ve gozlemcinin konum vektorlerini tanimlayalim

R = [R_boy*m.cos(omega),R_boy*m.sin(omega),0]
r = skaler(r_boy,r_bir)

Rr = topla(R,skaler(-1,r))
alfa = 90 - 180*m.acos(carp(Rr,r_bir)/boy(Rr))/m.pi

Dogu_iz = carp(Rr,fi_bir)          # Dogu yonune izdusum
Kuzey_iz = carp(Rr,teta_bir)      # Kuzey yonune izdusum
beta = 180*m.atan(Dogu_iz/Kuzey_iz)/m.pi

if Dogu_iz*Kuzey_iz >= 0:           # Bu kisim coklu deger alabilen
   if Dogu_iz < 0:                          # trigonometrik fonksiyonlarin
      beta = beta + 180                  # degerlerinin karismasini onlemek
else:                                              # icindir.
   if Kuzey_iz < 0:                          # 
      beta = beta + 180                  #
   else:                                           #
      beta = beta + 360                  #

print beta,'yonunde.'
print 'Ufuktan',alfa,'derece yuksekte.'

21 Mart 2012 Çarşamba

Yer-sabit uydular


Fizik yardımı ile ne kadar az bilgiden ne kadar çok şey üretilebileceğini gösteren şahane bir örnektir bu. Öyle ki sadece Newton yerçekimi kanununu bilen meraklı bir lise öğrencisi bu yazıdaki herşeyi tek başına üretebilir.

Dünya üzerindeki bir gözlemciye göre gökyüzünde SABİT bir noktada duran uydulara yer-sabit (geostationary) uydular denir. Televizyon yayınları başta olmak üzere birçok iletişim uygulamalarında aktif olarak kullanılmakta olan bu uyduların yörüngelerinin özelliklerinin ne olması gerektiği üzerinde düşünelim.

1. YÖRÜNGE EKVATOR DÜZLEMİNDE OLMALIDIR. Zira ekvator düzlemi ile 0'dan farklı bir açı yapan herhangi bir yörüngede dolanan uydunun gökyüzündeki konumu nereden bakılırsa bakılsın sürekli "yukarı" ve "aşağı" olarak değişecektir; sabit bir noktada durması mümkün değildir. Bu o kadar sağduyuya hitap eden birşey ki "matematiksel ispata" kalkışmak işleri kolaylaştıracağına zorlaştırır. :-))

2. YÖRÜNGE ÇEMBERSEL OLMALIDIR. Eliptik yörüngelerde dolanan uydular  bulundukları konuma göre hızlanır ve yavaşlarlar. (2. Keppler kanunu) Dünyanın dönme hızı ise (çok küçük değişimler ihmal edilirse) sabit olduğu için eliptik yörüngelerde olan uyduların yere göre sabit durması söz konusu olamaz.

3. UYDUNUN DÖNME PERİYODU DÜNYANINKİNE EŞİT OLMALIDIR. Elbette! Aksi takdirde ya geri kalır ya da ileri gider, iki durumda da gökyüzünde sabit duramaz.

Bu son koşul bize uydunun ne kadar uzakta olması gerektiğini de söyler çünkü periyot ile yarıçap arasında kesin bir bağıntı vardır. Bu bağıntı Newton yerçekimi kanununu kulanarak klasik mekanik ile kolayca hesaplanabilir. R yarıçaplı çembersel bir yörünge düşünelim. Dünya ile uydu arasındaki çekim kuvveti uyduyu böyle bir yörüngede tutmak için gerekli merkezcil kuvvete eşitlenirse:



Burada uydunun açısal hızı olup periyoda aşağıdaki şekilde bağlıdır:



Bu ifadeyi üstteki denklemde yerine yazıp uydunun kütlesini iki tarafta sadeleştirir ve R'yi yalnız bırakırsak aşağıdaki bağıntıyı elde ederiz.



Bu formülde dünyanın kütlesi (M), evrensel yerçekimi sabiti ve (G) dünyanın dönme periyodu olan 23 saat 56 dakika 4.1 saniye yerine yazılırsa yarıçap  42,164 km olarak bulunur.

Bu büyük çember üzerine yerleştirilen uyduların konumları ekvator üzerine izdüşümlerinin MERİDYENİ ile ölçülür. Şu anda yörüngede bulunan yer-sabit uyduların ve konumlarının listesine buradan ulaşılabilir. Görülebildiği gibi  uydular 0.1 derecelik farklarla bile yerleştirilebiliyorlar ki bu fark çember üzerinde yaklaşık 73 km. lik bir mesafeye denk gelir. Burada esas dikkat edilmesi gereken sinyallerin girişim yapıp yapmamasıdır ki bu teknik mevzuya girmiyoruz.

Bir sonraki yazıda vektör hesabı ve iki farklı koordinat sisteminin dönüşümlerini kullanarak dünyanın herhangi bir yerinden bakan bir gözlemcinin belli bir yer-sabit uyduyu hangi pozisyonda göreceğinin hesaplanmasını göstereceğim.

16 Mart 2012 Cuma

Nanotüp Nasıl İnşa Edilir?

Karbon nanotüplerin yapısı grafin (graphene) denilen iki boyutlu balpeteğine benzer karbon tabakalarının yukarıdaki şekilde gösterildiği gibi kıvrılıp bir silindir oluşturması şeklinde tasvir edilebilir. Hesaplamalı fizikte ve kimyada bu yapıların incelenmesi için öncelikle atomlarının koordinatlarının bilinmesi gereklidir. Geçen sene bir konferansda birisi bana bu işi yapabilen bir program bilip bilmediğimi sordu. Ben de hocamdan aldığım onun ABD'deki eski grubundaki birinin yazdığı programı fazla düşünmeden verdim. Ancak sonra böyle bir paylaşımın izinsiz yapılmasının sorun teşkil edebileceği söylendi bana. İnternette ve değişik kimya programlarında bu işi yapan hazır algoritmalar mevcut fakat açık kod olarak bulmak aktif olan bu konudaki rekabetten dolayı olsa gerek pek mümkün değil. Araştırmacılara faydalı olması ümidi ile burada bu iş için python dilinde yazdığım bir programın kodunu ve çalışma mantığını açıklıyorum.



Karbon nanotüpler "kiralite" (chirality) denen (m,n) gibi bir sayı çifti ile tasvir edilirler. Bu (m,n) sayı çifti grafinin baz vektörleri cinsinden kiralite vektörünü tanımlar. Örneğin yukarıdaki şekilde grafinin baz vektörleri a ve b ile gösterilmiş olup (4,2) kiralitesindeki bir nanotüp için kiralite vektörü (C) gösterilmiştir. a vektörünün 4 katı b vektörünün 2 katı ile toplandığı zaman C vektörünü verir.

Peki bu C vektörü ne işe yarar? C vektörü grafin tabakasının nasıl kıvrılıp tüp şekline getirilmesi gerektiğini gösterir. Tabaka öyle kıvrılmalıdır ki vektörün kuyruğundaki ve ucundaki örgü noktaları birbirinin üzerine gelsin. C haricinde önemli olan bir vektör de T vektörüdür. Bu vektör ismini ötelemeden (translation) alır ve tüpün ekseni boyunca kendini tekrarladığı minimum mesafeyi yani hücreyi tanımlar. Elbette ki T vektörü C vektörüne diktir. Nanotübümüzü inşa edebilmek için C'nin yanında T'yi de bilmek gerekir. O zaman işe T'yi hesaplamakla başlayalım. C'ye dik olma koşulunu yazarsak:



olmak zorundadır. Parantezleri açıp aralarında 60 derece olan a ve b'nin nokta çarpımlarını yaparsak:



olur ve bu da aşağıdaki şekilde yazılabilir:



Burada T_a ve T_b'nin en küçük değerleri bize lazımdır. O yüzden sağdaki kesri sadeleştirebildiğimiz kadar sadeleştirmemiz lazımdır. Bunu matematiksel olarak ifade etmek (ve bilgisayara anlatmak) için payın ve paydanın ortak bölenlerinin en büyüğünü (obeb) bulup bölmemiz lazımdır. Böylece T_a ve T_b için aşağıdaki eşitlikleri elde ederiz.




Şimdi buraya kadar olan kısmı bilgisayara anlatmaya çalışalım:

-------------------------------------------------------------
def obeb(A,B):
    while A:   
       A, B = B%A, A
    return B

m = input('Lutfen m\'yi giriniz: ')
n = input('Lutfen n\'yi giriniz: ')
hucre = input('Nanotup ekseni boyunca birim hucre sayisi: ')

T_a = -1*hucre*(2*n+m)/obeb(2*m+n,2*n+m)
T_b = hucre*(2*m+n)/obeb(2*m+n,2*n+m)
---------------------------------------------------------------

Elbette en küçük hücreli nanotüpü istemiyor olabiliriz, bu yüzden kullanıcıya kaç birim hücre istediğini de sorduk. En küçük birimi bulmak istiyorsak buraya 1 yazmamız lazım. Obeb için kullandığımız minik ve şirin fonksiyon ise Öklit algoritması olarak bilinir. T vektörünü hesapladığımıza göre devam edebiliriz. Kıvırma işlemi sonucunda koordinatlarını yazmamız gereken bölge yukarıdaki şekilde gri ile gösterilen C ve T'nin oluşturduğu dikdörtgendir. Bir sonraki aşama bu kısmın içerisinde yer alan örgü noktalarının k ve l indislerinin belirlenmesi olmalıdır. Mesela şekildeki K noktası (1,3) indislerine sahiptir (1a+3b). İşte meselemiz hangi indis çiftinin gri bölge içerisinde yer alacağıdır. Bunu yapmak için vektörel çarpıma başvuruyoruz. Seçilen herhangi bir K noktasının içeride kalması için aşağıdaki eşitsizliklerin hepsinin sağlanması gerekmektedir:






Böyle birşey yazdığım için matematikçiler bana kızabilirler zira vektörlerle sayılar kıyaslanamaz ama burada kastettiğim şey sonuç vektörünün sayfanın dışına doğru mu içine doğru mu baktığıdır. Yani z-koordinatlarının işaretleridir. Parantezleri açıp vektörel çarpımları yaptığımızda bu ifadeyi aşağıdaki iki eşitsizlik haline getirebiliriz:




Şimdi bu söylediklerimizi koda dökme zamanı. Aklımda olan şey k'ya ve l'ye belli aralıktaki değerler için döngü yaptırıp yukarıdaki koşulu sağlayan çiftleri bir listeye kaydetmek. Tabii aralık ne olacak sorusu var; bunun için m,n,T_a ve T_b indislerinin toplamı yeterli olması lazım.
-----------------------------------------------------

aralik = m + n + abs(T_a) + abs(T_b) + 2
kl_listesi = [ ]

for k in range(-aralik,aralik):
    for l in range(-aralik,aralik):
        if (m*l-n*k) >= 0 and (m*T_b-n*T_a) > (m*l-n*k):
           if (k*T_b-T_a*l) >= 0 and (m*T_b-n*T_a) > (k*T_b-T_a*l):
              kl = [ ]
              kl.append(k)
              kl.append(l)
              kl_listesi.append(kl)

------------------------------------------------------

Böylece kl_listesi'ni ihtiyacımız olan ikililerle doldurmuş oluruz. Şimdi artık atomları yerleştirmeye (koordinatlarını yazmaya) başlayabiliriz. Bunun için elbette birim hücrede yer alan atomları önceden bir listeye yazmak gerek. Grafin'in birim hücresinde iki karbon atomu bulunur. Bunlar şekildeki içi dolu ve boş noktalarla gösterilmiştir. Yine içiçe döngü ile her baz atomunu kl_listesindeki her ikili ile eşleştirip ötelemek lazımdır. Tabii bir de baz vektörlerini de artık yazmanın vakti geldi. Bunu da grafindeki karbon-karbon bağı uzunluğu olan 1.42 Angstrom cinsinden yazabiliriz.

-------------------------------------------------------
bag_uzunlugu = 1.42

a_x = bag_uzunlugu*(3**0.5)
a_y = 0.0
b_x = a_x/2
b_y = b_x*(3**0.5)

baz_atomlari = [[0,0],[b_x,bag_uzunlugu/2]]

atom_listesi = [ ] 
for i in kl_listesi:
    for j in baz_atomlari:
        koordinatlar = [ ]
        koordinatlar.append(j[0]+i[0]*a_x+i[1]*b_x)
        koordinatlar.append(j[1]+i[0]*a_y+i[1]*b_y)
        atom_listesi.append(koordinatlar)

-----------------------------------------------------

Artık gri dikdörtgendeki atomların x-y koordinatları atom_listesi'nin içinde. Şimdi bu dikdörtgeni kıvırıp bir nanotüp yapmak istiyoruz. Bu işi yapan programlar genelde nanotüp eksenini z-ekseni ile çakışık hale getirirler. Bunu yapmanın en kolay yolu x-y eksenine göre teta açısı kadar yan duran dikdörtgenimizi "düzeltmek"ten geçer. Düzeltme ile kast ettiğimiz şey bir koordinat dönüşümüdür. Bir koordinat sistemini orjini etrafında saat yönünün tersinde teta açısı kadar döndürdüğünüzde yeni koordinatlar aşağıdaki formülle verilir.



Biz koordinat sistemimizi saat yönünde döndürmek istiyoruz dolayısı ile teta'nın işareti negatiftir bu da sadece sinüslerin işaretini ters çevirir. Elbette ki önce teta açısını hesaplamamız gerekir, bunun için de C vektörünün bileşenlerine ihtiyacımız vardır. Bir de trigonometrik fonksiyonları kullanma için python'un math modülünü çağırmamız gerekmektedir.

--------------------------------------------------

C_x = a_x*m + b_x*n
C_y = a_y*m + b_y*n

C_boy = (C_x**2 + C_y**2)**0.5

import math

teta = math.acos(C_x/C_boy)

dondurulmus = [ ]
for i in atom_listesi:
    koordinatlar = [ ]
    koordinatlar.append(i[0]*math.cos(teta)+i[1]*math.sin(teta))
    koordinatlar.append(-1*i[0]*math.sin(teta)+i[1]*math.cos(teta))
    dondurulmus.append(koordinatlar)

---------------------------------------------------



Artık dikdörtgenimizin içinde yer alan atomlar yukarıdaki şekildeki gibi kayıt edilmişlerdir. (Dikey eksen henüz halen y-eksenidir ama bir sonraki basamakta y koordinatını z-koordinatı olarak yazacağız, burada yandaki şekille karışmasın diye erkenden z olarak isimlendirdik.) Son basamak olan kıvırma işlemini yapabilmek için bütün atomları yanda gösterildiği gibi bir çemberin üzerine dizmek gereklidir. C vektörünün uzunluğu çemberin çevresidir ve A'nın x koordinatının C'nin boyuna oranının 2 pi ile çarpımı alfa açısını radyan cinsinden bize verir. x ve y koordinatları alfanın trigonometrik oranları cinsinden yazılabilir:

------------------------------------------------------
T_x = a_x*T_a + b_x*T_b
T_y = a_y*T_a + b_y*T_b
T_boy = (T_x**2 + T_y**2)**0.5

R = C_boy/(2*math.pi)

kivrilmis = [ ]

for i in dondurulmus:
    koordinatlar = [ ]
    alfa = 2*math.pi*i[0]/C_boy
    koordinatlar.append(R*math.cos(alfa))
    koordinatlar.append(R*math.sin(alfa))
    koordinatlar.append(i[1]-T_boy/2)  # orjini z boyunca ortalamak icin
    kivrilmis.append(koordinatlar)

---------------------------------------------------------

Artık sadece yazdırmak kaldı. Bunu birçok kimya programının tanıdığı standart .xyz dosyası formatında yapalım. Bu formatta önce başa toplam atom sayısı yazılır, sonraki satır açıklamalara ayrılmıştır (boş bırakılabilir) ve sırasıyla atom sembolü, x,y ve z koordinatları araya boşluk bırakarak yazılır.

---------------------------------------------------------

print len(kivrilmis)
print
for i in kivrilmis:
    print '%s  %8f  %8f  %8f'  % ('C',i[0],i[1],i[2])

----------------------------------------------------------

Çıktıyı bir dosyaya yazdırıp MOLEKEL programı ile çizdirdiğim 4 birim hücreli bir (5,5) nanotüpün resmi aşağıdadır.



Evet! Bu kısa kodla ilk resimde gösterilen işlemi yapmış olduk. Aynı mantıkla çalışan bir kodla prensipte sadece grafinden yuvarlanarak yapılan nanotüpler değil her türlü peryodik düzlemsel yapıdan elde edilebilecek nanotüplerin koordinatlarını bulmak mümkündür. Aynı "kiralite" mantığı birçok farklı yapıya genelleştirilip kullanılabilir.

Yapılması gereken değişiklik baz vektörlerinde ve baz atomlarındadır. Bir de T vektörünü indisler cinsinden bulurken kodu yazmadan önce matematiği yukarıda gösterdiğimiz metod ile elle yapmak ve tam sayılar cinsinden bir formül türetmeye çalışıp ondan sonra koda yazılmalıdır zira bu işi baz vektörleri bileşenlerinden bilgisayarlara yaptırmaya kalkmak "floating point" meseleleri yüzünden riskli ve zahmetlidir. Yine gerekli indisleri yazdırırken kullanılan vektörel çarpımlar yeni baz vektörlerine göre kontrol edilmelidir. En kritik hususlar bunlardır. 

Bunlar haricinde nanotüpü eksen boyunca germe, sıkıştırma, eksene dik ölçeklendirme (şişirme veya büzüştürme) veya eksen etrafında burma küçük modifikasyonlarla elde edilebilir. Bu modifikasyonların da eklendiği en genel programı aşağıda veriyorum. ### ile işaretli yerler grafin dışında bir yapı ile çalışılmak isteniyorsa değiştirilmesi gereken yerlerdir. Son olarak seçtiğiniz baz vektörlerinin axb vektörel çarpımının pozitif olması lazımdır.

-----------------------------------------------------------------
import math   # Trigononetrik fonksiyonlar icin lazim

def obeb(A,B):    # Ortak bolenlerin en buyugu fonksiyonu
    while A:
       A, B = B%A, A
    return B

### bag_uzunlugu = 1.42

### a_x = bag_uzunlugu*(3**0.5)
### a_y = 0.0
### b_x = a_x/2
### b_y = b_x*(3**0.5)


### baz_atomlari = [['C',0,0],['C',b_x,bag_uzunlugu/2]]

print 'Kiralite (m,n) olarak tanimlaniyor'
m = input('m\'yi giriniz: ')
n = input('n\'yi giriniz: ')
hucre = input('Hucre sayisini giriniz: ')
T_olcek = input('T olceklendirmesi: ')
C_olcek = input('C olceklendirmesi: ')
burulma = input('Burulma acisi (derece): ')

### T_a = -1*hucre*(2*n+m)/obeb(2*m+n,2*n+m)
### T_b = hucre*(2*m+n)/obeb(2*m+n,2*n+m)

aralik = m + n + abs(T_a) + abs(T_b) + 2
kl_listesi = [ ]

for k in range(-aralik,aralik):
    for l in range(-aralik,aralik):
        if (m*l-n*k) >= 0 and (m*T_b-n*T_a) > (m*l-n*k):
           if (k*T_b-T_a*l) >= 0 and (m*T_b-n*T_a) > (k*T_b-T_a*l):
              kl = [ ]
              kl.append(k)
              kl.append(l)
              kl_listesi.append(kl)

atom_listesi = [ ]
for i in kl_listesi:
    for j in baz_atomlari:
        koordinatlar = [ ]
        koordinatlar.append(j[0])
        koordinatlar.append(j[1]+i[0]*a_x+i[1]*b_x)
        koordinatlar.append(j[2]+i[0]*a_y+i[1]*b_y)
        atom_listesi.append(koordinatlar)

C_x = a_x*m + b_x*n
C_y = a_y*m + b_y*n
C_boy = ((C_x**2 + C_y**2)**0.5)

teta = math.acos(C_x/C_boy)
dondurulmus = [ ]

for i in atom_listesi:
    koordinatlar = [ ]
    koordinatlar.append(i[0])
    koordinatlar.append(i[1]*math.cos(teta)+i[2]*math.sin(teta))
    koordinatlar.append(-1*i[1]*math.sin(teta)+i[2]*math.cos(teta))
    dondurulmus.append(koordinatlar)

T_x = a_x*T_a + b_x*T_b
T_y = a_y*T_a + b_y*T_b
T_boy = T_olcek*((T_x**2 + T_y**2)**0.5)

R = C_olcek*C_boy/(2*math.pi)
kivrilmis = [ ]

for i in dondurulmus:
    koordinatlar = [ ]
    alfa = 2*math.pi*i[1]/C_boy
    beta = (T_olcek*i[2]*burulma/T_boy)*(math.pi/180)
    koordinatlar.append(i[0])
    koordinatlar.append(R*math.cos(alfa+beta))
    koordinatlar.append(R*math.sin(alfa+beta))
    koordinatlar.append(T_olcek*i[2]-T_boy/2)
    kivrilmis.append(koordinatlar)

cikti = [ ]
for i in kivrilmis:
    a = '%s  %8f  %8f  %8f\n' % (i[0],i[1],i[2],i[3])
    cikti.append(a)

isim = raw_input('Lutfen cikti dosyasi icin bir isim giriniz: ')
dosya = file(isim,'w')
dosya.write(str(len(atom_listesi)))
dosya.write('\n')
a = '('+str(m)+','+str(n)+') Nanotup. R = '+str(R)+' Ang. T = '+str(T_boy)+' Ang.\n'
dosya.write(a)

for i in cikti:
    dosya.write(i)

7 Mart 2012 Çarşamba

Rastgelenin gücü: Monte Carlo



Pi sayısını hesaplamanın değişik ve sıradışı bir yöntemi vardır:

İçerisine çeyrek çember çizilmiş bir karenin içinde RASTGELE noktalar seçelim. Eğer seçimleriniz gerçekten rastgele(!) ise çeyrek çemberin içinde kalan noktaların sayısının tüm noktaların sayısına oranının çeyrek çemberin alanının karenin alanına oranına yakın bir sayı olmasını bekleriz. Bu alanların oranı da kolayca görülebileceği gibi 'e eşittir. Dolayısı ile sadece bir rastgele sayı üreteci ve bir koşul belirleyerek insanlığın bu en eski problemlerinden biri nispeten "zahmetsizce" çözülebilir. 

Bu tipteki yaklaşımlara Monte Carlo Yöntemi ismi verilir. 20. yy'nın büyük hezarfeni John von Neumann tarafından Los Alamos'da atom bombası üzerinde çalışırken keşfedildiği rivayet edilir ve olasılığa (kabaca "şansa") dayalı yapısı nedeni ile meşhur kumarhanenin ismi ile isimlendirilmiştir. 

Günümüzde Monte Carlo yöntemi fizikten mühendisliğe, hesaplamalı biyolojiden uygulamalı istatistiğe, oyunlardan dizayna ve görsel uygulamalara, finans ve iş dünyası modellerinden telekominikasyona kadar geniş bir yelpazede kendisine uygulama alanı bulur. 

Sadece yukarıdaki pi sayısı örneğini düşünürsek üstteki paragraf biraz abartılı gelebilir zira bu problem az değişkenli bir problemdir ve geleneksel yollardan rahatça çözülebilir. Monte Carlo'nun gücü kendisini çok değişkenli fonksiyonlarda gösterir çünkü bu fonskiyonların integral hesaplamalarının klasik yollarla alacağı zaman fonksiyonun değişken sayısına ÜSTEL bağımlılık gösterir.  Oysa Monte Carlo algoritması "boyutsallık laneti" diye isimlendirilen bu durumdan bu kadar ağır etkilenmez ve yüksek boyutlu problemlerde avantaj sağlamaya başlar.

Yukarıda animasyonu verilen pi sayısı hesaplama problemi için python dilinde yazdığım kısacık bir bilgisayar kodu ile bitirelim. Deneyip kendiniz de görebilirsiniz.

import random        # Python'da rastgele sayi üreten protokolun yer aldigi modul
nokta_sayisi = 100000      # Toplam nokta sayisi
icinde = 0                          # cemberin icinde kalan nokta sayisi
for m in range(nokta_sayisi):
    x = random.random()       # bir noktanin x koordinati
    y = random.random()       # ayni noktanin y koordinati
    if y < (1-x**2)**0.5:    # Cember denkleminden cemberin icinde kalma kosulu
       icinde += 1          # eger kosul dogruysa icindeki nokta sayisini 1 arttirir
pi_sayisi = 4*float(icinde)/float(nokta_sayisi)
print 'Pi sayisi yaklasik olarak:', pi_sayisi

Programı yürüttüğümde aldığım bazi sonuclar:
3.14096
3.14396
3.14684
3.13652
3.1386
3.13224  vs. vs. vs....

3 Mart 2012 Cumartesi

Kainatın en temel ilkesi


Kendimizi üstteki şeildeki cankurtaranın yerine koyalım. Denizde bir adam boğulduğunu görüyorsunuz ve yerinizden fırlayıp adamı kurtarmaya koşuyorsunuz. Hangi yolu takip edersiniz? EN KISA olanı mı? Yoksa sizi adama EN KISA ZAMANDA ulaştıracak olanı mı? Elbette ikincisi. Koşma hızınız yüzme hızınızdan fazla olduğuna göre  bunu yapmak için hangi noktaya kadar koşup ondan sonra yüzmeye başlamamız gerektiği üzerinde düşünelim. Suya girdiğimiz nokta O olsun. Karadaki hızım sudaki hızım ise karada geçirdiğim zaman iken suda geçirdiğimiz zaman olur.Bu zaman ifadelerini x, h1, h2 ve L cinsinden de yazıp toplarsak toplam varış zamanımızı aşağıdaki gibi yazabiliriz:



İşte bu zamanı en kısa yapacak x'i bulmak istiyorum. Yani problemimiz t(x) fonksiyonunun minimumunu bulma problemine indirgenmiş oldu, bunu da x'e göre türev alıp sıfıra eşitleyerek yapabiliriz. Bu türevi alıp sıfıra eşitlersek aşağıdaki denklemi bulmuş oluruz.



Bu denklemi çözersek prensipte x'i bulabiliriz ancak buna kalkışmayacağım zira burada vurgulamak istediğim başka bir şey. v1 ve v2'yi ayrı tutarsak geride kalan terimlerin şekilde (ne için şaretlediğimi merak ettiğiniz) açıların sinüsleri olarak yazılabileceği temel dik üçgen bağıntılarından görülebilir. Denklemimizi yeniden düzenlersek aşağıdaki şekle getirebiliriz.



Tanıdık geldi mi biryerlerden?



Evet, ışığın kırınımındaki Snell yasasından bahsediyorum:



Buradaki n kırınım indisi denilen ve ışığın ortamdaki hızı ile ters orantılı olan, ortama bağlı bir sabittir. Görüldüğü gibi Snell yasası üstte türettiğimiz formülden başka birşey değil. Yani Maxwell denklemlerini sınırda çözmek yerine IŞIK EN KISA ZAMANDA VARACAK ŞEKİLDE bir yol izler prensibini kabul edersek Snell yasasını çok daha "kısa, basit ve değişik" bir yoldan türetebiliriz. Bu prensibe Fermat prensibi denir.

"Güzel bir oyuncak ama başlık biraz iddialı değil mi? Sonuçta sadece ışığın hareketi hakkında bir ilke bu" diyebilirsiniz ancak Fermat prensibi esasında çok daha genel bir prensibin özel bir halidir. Bahsettiğimiz bu ilke EN AZ EYLEM ilkesidir. (Principle of least action). Bu ilke, kainattaki değişimlerin her zaman "eylem"i en az yapacak şekilde gerçekleşeceğini söyler. (Ben bunu kainat her zaman en "kestirme yolu" izler diye ifade etmeyi seviyorum.) Burada özel bir anlam yüklediğimiz "eylem" ile ne kastettiğimizi bir örnek ile açalım. Klasik mekanik üzerinden konuşacak olursak eylem (S) sistemin kinetik enerjisi (T) ile potansiyel enerjisinin (U) farkının zaman integrali olarak tanımlanır.



Yani t1 ve t2 zaman aralığında sistem öyle bir değişir ki yukarıdaki fonksiyon en küçük değeri alsın. Bu kadar! Ne eylemsizlik prensibi, ne F = ma, ne etkiye tepki hiçbir Newton prensibini kullanmadan  sadece azıcık matematik bilerek tamamıyla eşdeğer (ve esasında daha genel) bir klasik mekanik kurmak mümkündür! İnanılmaz değil mi!

Ünlü fizikçi Feynman bu mevzu hakkında şöyle söylüyor: "Lise hocam bir gün fizik dersinden sonra beni yanına çağırdı ve "derste seni sıkılmış görüyorum; gerçekten ilginç birşey duymak ister misin" dedi ve öyle bir şey söyledi ki bunu inanılmaz derecede büyüleyici buldum. O günden beridir de hala büyüleyici bulurum. Söylediği şey "en az eylem" ilkesiydi."

En az eylem ilkesinin uygulaması sadece klasik mekanik ile de sınırlı değildir. Klasik alan teorileri, kuantum mekaniği, kuantum alan teorisi, izafiyet teorisi gibi çok geniş sahalara genelleştirilmiştir ve bu, fiziksel bilimlerin en önemli genelleştirmelerinden biridir.

Evet, kainat her zaman en kestirme yolu izler...